LATE QUATERNARY GLACIER FLUCTUATIONS,
ITSI RANGE, YUKON TERRITORY
by
Caleb Andrew Mathias
B.Sc., University of Tennessee, 2019

THESIS SUBMITTED IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF
MASTER OF SCIENCE
IN
NATURAL RESOURCES AND ENVIRONMENTAL STUDIES

UNIVERSITY OF NORTHERN BRITISH COLUMBIA
April 2025
© Caleb Andrew Mathias, 2025

Abstract
Understanding past glacier fluctuations is critical for understanding Holocene
paleoclimate, yet glacier chronologies in the northeastern Canadian Cordillera remain sparse and
are frequently obscured by the occurrence of geologic scatter. This thesis reconstructs late
Pleistocene and Holocene glacier changes in the Itsi Range, Yukon Territory, and examines
processes driving regionally consistent scatter in moraine ages. Using terrestrial in situ 10Be
surface exposure dating, I identify two Little Ice Age advances at the Itsi Glacier among the Itsi
Range: 0.96 ka ± 0.36 (outer moraine) and 0.20 ka ± 0.05 (inner moraine). Erratics dated to
11.61 ka ± 0.54 indicate Cordilleran Ice Sheet retreat concurrent with the Younger Dryas
termination. To contextualize the Itsi Glacier chronology with other regional advances, I
compiled 10Be dated moraine records across western North America moraines organized by
moraine morphostratigraphic position as denoted by end groups 1 – 3. The Itsi Glacier inner
moraine aligns closely with end group 1 (0.13 ka ± 0.10), while outer moraine ages reflect
significant scatter, mirroring the poorly constrained end group 2 record (0.46 ka ± 0.35).
Distinctions among end group 1 and 2 moraine ages across latitudes suggest that northeastern
Canadian Cordillera glaciers emplaced their moraines centuries earlier than those in the
southwest. I assess the western North America moraine record with relevant multicentennialscale climate proxy records to deduce what factors influence glacier growth and retreat. Late
Holocene advances were driven by orbital forcing, with additional contributions from
temperature, precipitation, greenhouse gas emissions, and solar irradiance anomalies coherent
with end group 1 culminations. End group 2 moraine scatter complicates its climatic
interpretation. A parsimonious explanation for the inconsistencies among end group 2 moraine
ages point to the result of amalgamated moraines, heterogeneous mixtures of previously exposed

ii

and reincorporated moraine boulders, deposited between 3.54 – 2.50 ka during the TiedemannPeyto Advance. Further research must be carried out to quantify geomorphic uncertainties
inherent to techniques of cosmogenic nuclide dating.

iii

Table of Contents
Abstract ............................................................................................................................... ii
Table of Contents ............................................................................................................... iv
List of Tables .................................................................................................................... vii
List of Figures .................................................................................................................. viii
Acknowledgements ............................................................................................................ xi
1.

2.

Introduction ................................................................................................................. 1
1.1

Motivations and objectives ................................................................................. 2

1.2

Document structure ............................................................................................. 2

Holocene activity among the Itsi Range ..................................................................... 3
2.1

Climate and the cryosphere ................................................................................. 3

2.1.1

Spatial and temporal climate scales ................................................................ 3

2.1.2

Glacier chronology as climate proxy .............................................................. 5

2.1.3

Modern cryospheric conditions....................................................................... 7

2.1.4

Canadian Cordillera cryospheric activity........................................................ 8

2.2

Terrestrial in situ cosmogenic nuclide dating theory ........................................ 10

2.3

Study site........................................................................................................... 13

2.4

Methodology ..................................................................................................... 17

2.4.1

10

2.4.2

DEM creation and geodetic mass balance .................................................... 21

Be dating .................................................................................................... 17

iv

2.4.3
2.5

Termini tracking............................................................................................ 26
Results ............................................................................................................... 28

2.5.1

Itsi Glacier 10Be chronology ......................................................................... 28

2.5.2

Termini tracking............................................................................................ 33

2.5.3

Geodetic mass balance .................................................................................. 33

2.6

Discussion ......................................................................................................... 34

2.6.1

Geologic scatter ............................................................................................ 34

2.6.2

Inner moraine ................................................................................................ 39

2.6.3

Outer moraine ............................................................................................... 40

2.6.4

Erratics .......................................................................................................... 41

2.6.5

Geodetic mass balance and termini retreat ................................................... 42

2.7
3.

Conclusion ........................................................................................................ 46

Comparison with western North America moraine records ..................................... 47
3.1

Introduction ....................................................................................................... 47

3.2

Methodology ..................................................................................................... 48

3.3

Discussion of moraine records .......................................................................... 51

3.3.1

Review of western North America Neoglacial advances ............................. 51

3.3.2

Radiocarbon, dendrochronology, and glaciolacustrine records .................... 53

3.3.3

10

3.3.4

10

Be moraine ages ......................................................................................... 55
Be erratic ages ............................................................................................ 61
v

4.

5.

3.4

Limitations ........................................................................................................ 62

3.5

Conclusion ........................................................................................................ 65

Comparison of moraine record with climate proxy records ..................................... 67
4.1

Introduction ....................................................................................................... 67

4.2

Methodology ..................................................................................................... 68

4.3

Results and discussion ...................................................................................... 73

4.3.1

Moraine end groups ...................................................................................... 74

4.3.2

Anthropogenic warming ............................................................................... 76

4.4

Limitations ........................................................................................................ 77

4.5

Conclusion ........................................................................................................ 79

Conclusions ............................................................................................................... 80
5.1

Summary of advances ....................................................................................... 80

5.2

Research limitations .......................................................................................... 82

5.3

Research recommendations .............................................................................. 83

Appendix ........................................................................................................................... 85
References ......................................................................................................................... 91

vi

List of Tables
Table 2.1: Itsi Glacier field data collected on 13 August 2022. ................................................... 18
Table 2.2: Aerial photographs used for glacier change analysis. ................................................. 22
Table 2.3: NMAD errors before and after coregistration for the DEM of each year. .................. 24
Table 2.4: Itsi Glacier 10Be sample information. .......................................................................... 32
Table 2.5: Statistics derived from the GTT tool and area analysis among all glaciers of the Itsi
Range. ........................................................................................................................................... 33
Table 2.6: Mass budget of Itsi Range glaciers (1949 – 2022). ..................................................... 34
Table 2.7: Itsi Glacier moraine and erratic ages with accompanying reduced chi-squared (X2)
and Mean Absolute Deviation (MAD) statistics. .......................................................................... 36
Table 3.1: Description for each morphostratigraphic position used for WNA comparisons. ....... 49
Table 3.2: 10Be dated WNA moraines grouped by morphostratigraphic position with associated
statistics (Dortch et al., 2010; Hawkins et al., 2023, 2021; Howley, 2008; Menounos et al.,
2013a; Pendleton et al., 2017; Young et al., 2009). ...................................................................... 57

vii

List of Figures
Figure 2.1: July solar insolation at 65° N (W/m2) (Berger and Loutre, 1991) and NH 50 year
ensemble mean temperature anomalies (Erb et al., 2022) throughout the Holocene. .................... 4
Figure 2.2: (a) Regional view of the study site in relation to adjacent physiographic provinces.
(b) Itsi Range study area and glaciers overlayed on Sentinel 2 10 m true colour imagery captured
on 11 August 2022. ....................................................................................................................... 17
Figure 2.3: Three-dimensional visualization of the 2004 orthomosaic of the Itsi Range, optimized
using CLAHE. No vertical exaggeration is applied. Glaciers are outlined in red, outer moraine
highlighted in yellow, and inner moraine highlighted in blue. Scale is approximate. Created in
QGIS using the Qgis2threejs plugin. ............................................................................................ 23
Figure 2.4: Itsi Glacier and moraines. Aerial photo and Sentinel-2 derived glacier extents are
overlayed for 1949, 2004, and 2022. DEM basemap was generated using stereoscopic aerial
image pairs from 2004. Individually 10Be dated boulder ages from each moraine are included and
reported with propagated error as LSDn exposure age BP. .......................................................... 29
Figure 2.5: 10Be ages from two neoglacial moraines emplaced by the Itsi Glacier reported in ka
BP. Grey curves show ages of individual boulders rendered as Gaussian distributions. The
coloured lines show the sum of each individual Gaussian distribution for the 10Be ages sorted by
morphostratigraphic position. Box-and-whisker plots corresponding to each moraine denote the
median ± IQR, minimum, and maximum surface exposure ages. ................................................ 30
Figure 2.6: 10Be ages from two erratic boulders outboard the Itsi Glacier moraines reported in ka
BP. The pink line shows the sum of each individual Gaussian distribution. Box-and-whisker plot
denotes the median ± IQR, minimum, and maximum surface exposure ages of the two erratics. 31
Figure 2.7: A regional overview of permafrost zone coverage in relation to the Itsi Range field
site from Obu et al. (2019). ........................................................................................................... 38
Figure 2.8: An oblique angle of the Itsi Glacier forefield highlighting the inner and outer
moraines as well as an elongated ridge of glacially abraded rock proximal to the ice margin and
parallel to the moraine axes. ......................................................................................................... 41
Figure 2.9: ERA5-Land values for the Itsi Range (62.9°, -130.2°) plotted to represent boreal
summer melt (June, July, August (JJA)) and winter accumulation (December, January, February
(DJF)) seasonal means. Values are plotted as five-year moving averages with associated
trendlines. ...................................................................................................................................... 43
Figure 2.10: ERA5-Land 2023 MAT anomalies for the YT using reference period 1950-1980.
Itsi Range is highlighted in green. ................................................................................................ 44
Figure 2.11: Observed mean discharge values at Ross River (61.99, -132.41; Station ID:
09BA001) within the Upper Pelly sub-sub-drainage area for May and June. Extracted from the
Environment and Climate Change Canada Historical Hydrometric Data web site
(https://wateroffice.ec.gc.ca/mainmenu/historical_data_index_e.html). ...................................... 45
Figure 3.1: Regional compilation of glacial advance intervals from 4.0 ka to present based on
regional reviews using evidence from radiocarbon dating, dendrochronology, and
glaciolacustrine sedimentation records (Barclay et al., 2009; Hawkins et al., 2021; Menounos et
viii

al., 2009; Mood and Smith, 2015; Tomkins et al., 2008). Grey box widths are set by duration of
advance interval with ages calibrated to ka BP 1950. Dashed lines indicate approximate age
intervals when discrete ages were not provided. Orange and blue vertical bars are Itsi Glacier
median moraine ages, signifying moraine stabilization and glacier retreat; bar widths are set by
the IQR of 10Be dated moraine ages. Pink vertical bars denote ages of individual Itsi Glacier
outer moraine boulders which comprise the median outer moraine age; bar widths are set by the
propagated error of each boulder. ................................................................................................. 53
Figure 3.2: Regional compilation of WNA late Holocene 10Be boulder ages, including those
collected for this thesis, sampled from moraines and grouped by morphostratigraphic position
(Dortch et al., 2010; Hawkins et al., 2023, 2021; Howley, 2008; Menounos et al., 2013a;
Pendleton et al., 2017; Young et al., 2009). Thin grey curves show ages of individual boulders
rendered as Gaussian distributions. Thick curves represent the sum of each Gaussian distribution
for the entire population of 10Be moraine group ages, organised by colour. Itsi Glacier moraine
groups are inclusive in the regional curves and also highlighted separately within the summed
Gaussian distributions. Corresponding box-and-whisker plots denote the median ± IQR,
minimum, and maximum surface exposure ages of each moraine group. .................................... 56
Figure 3.3: Regional compilation of WNA late Holocene 10Be moraine ages, including those
collected for this thesis, grouped by morphostratigraphic position (Dortch et al., 2010; Hawkins
et al., 2023, 2021; Howley, 2008; Menounos et al., 2013a; Pendleton et al., 2017; Young et al.,
2009). Left: map illustrating location of sampled moraines where size of marker is proportional
to IQR. Right: scatter plot showing median moraine ages in relation to IQR. ............................. 59
Figure 3.4: Conceptual model of moraine amalgamation adapted from Bierman et al. (2021). An
amalgamated moraine is an unstratified mixture of clasts, some of which have experienced
periods of glacier transport and deposition one or several times before final transport and
deposition on the contemporary moraine. At time step 𝑡1, (i) buried clasts are shielded from
cosmic rays, (ii) glacially transported, and (iii) emplaced as a terminal moraine following
deglaciation. At 𝑡2, (iv) buried boulders are, (v) glacially transported, and (vi) emplaced as a
recessional moraine proximal to the ice margin. At 𝑡3, the glacier (vii) overrides the prior
moraine deposit from the second glaciation and (viii) emplaces an unstratified mixture of
previously deposited clasts as an amalgamated terminal moraine. The resulting exposure age
distribution (ix) reflects multiple records of deglaciation............................................................. 61
Figure 3.5: Regional compilation of WNA 10Be erratic boulder ages, including those collected
for this thesis, from BC and YT (Hawkins et al., 2023; Menounos et al., 2017; Stroeven et al.,
2010). Thin grey curves show ages of individual boulders rendered as Gaussian distributions.
The thick pink curves represent the sum of each Gaussian distribution for the entire population of
10
Be erratic ages. Itsi Glacier erratics are highlighted separately within the summed Gaussian
distributions. A corresponding box-and-whisker plot denotes the median ± IQR, minimum, and
maximum surface exposure ages for all erratics. .......................................................................... 62
Figure 4.1: Neoglacial climate forcings and indicators. (a) Global volcanic eruptions per decade
(30-year moving average) and major VSSI event magnitudes (Tg S) (Sigl et al., 2021; Toohey
and Sigl, 2019). (b) NPI at Mt. Logan (750-year spline) (Osterberg et al., 2014). (c) TSI
anomalies (ΔW m-2) plotted as 50-year moving averages (light line) and smoothed using a
Gaussian filter (dark line) (Steinhilber et al., 2009). (d) A spline smoothed dataset of atmospheric
GHGs CO2 and CH4 plotted as radiative forcing (ΔW m-2) (Köhler et al., 2017) Increases in
ix

GHG concentrations past 0.1 ka are not plotted. (e) July mean Solar Insolation at 65° N (W m-2)
at a 1000 year resolution (Berger and Loutre, 1991). (f) Itsi Range temperature (ΔT) and
precipitation anomalies (ΔP) (Liu et al., 2009) plotted as 50-year moving averages. (g) Regional
compilation of late Holocene 10Be boulder ages, including those collected for this thesis (Dortch
et al., 2010; Hawkins et al., 2023, 2021; Howley, 2008). Blue and orange lines represent the sum
of each 2σ Gaussian distribution for end groups 1 and 2, respectively. Itsi Glacier moraine end
groups are highlighted separately within the distributions. .......................................................... 73
Figure 4.2: Projected Itsi Range (62.9°N, -130.2°W,) SSP climate change scenarios (Cannon et
al., 2015; ClimateData.ca, 2024; McKenney et al., 2011; Tam et al., 2019)................................ 77

x

Acknowledgements
I acknowledge that the University of Northern British Columbia resides on the unceded
territory of the Lheidli T'enneh, the First Nation where the two rivers flow together. I thank the
Yukon Kaska Dena First Nations of the Ross River Dena Council and Liard First Nation for
allowing the sampling of boulders on their land. Thank you to Gerald Osborne and Rachel
Pagdin for accompanying Brian and me during field sampling.
I extend my deepest gratitude to my supervisor, Brian Menounos, for providing not only
the opportunity for this research but also what became the first vertiginous steps toward a new
life abroad. You offered invaluable guidance and provided opportunities that enriched my
research skills throughout this journey. Thanks to Adam Hawkins for his insights as both a friend
and collaborator during the development of my research, and to Charlie Bourque for keeping me
accountable through to the end. Lastly, I thank my partner, Parker Ahlstrom, for providing a
warm home and daily support, which were essential to completing my work.
This research was made possible by funding from NSERC and the Hakai Institute.

xi

1. Introduction
The principal goal of this thesis is to reconstruct paleoclimate behaviour during the
Holocene in the northeastern Canadian Cordillera, specifically among the Selwyn Mountains of
the Yukon Territory (YT). In this region, glacier chronologies remain sparse due to the remote
nature of the region and a lack of available organic materials for radiocarbon dating.
The timing of mountain glacier advances in the northeastern Canadian Cordillera, which
encompasses mountain ranges in the eastern YT and western Northwest Territories (NWT), in
response to climate events throughout the Holocene remains poorly understood (Hawkins et al.,
2021; Menounos et al., 2009). Refining glacier chronologies in this area will allow for improved
understanding of regional Holocene climate drivers and the complex morphological processes
acting upon moraines (Goehring et al., 2022; Solomina et al., 2015). A widely utilized method
for developing a glacier chronology is surface exposure dating through the measurement of
Terrestrial in situ Cosmogenic Nuclides (TCNs) (Balco, 2020a, 2011; Gosse and Phillips, 2001;
Schaefer et al., 2022; Solomina et al., 2015). TCN Dating (TCND) is a method to determine how
long a rock has been exposed to the flux of cosmic rays emitted primarily from within the Milky
Way galaxy (Gosse and Phillips, 2001). Certain minerals react to these cosmic rays with the
production of rare radioactive isotopes, the concentration of which is proportional to surface
exposure duration (Balco, 2011; Darvill, 2013; Gosse and Phillips, 2001)
Few applications of TCND have been carried out resolving Little Ice Age (LIA)
Holocene mountain glacier fluctuations in the northeastern Canadian Cordillera (Menounos et
al., 2009; Solomina et al., 2015). As glacier-fed catchments in this region face declines in latesummer streamflow and heightened risks to geomorphic stability and water quality,

1

understanding the response of alpine glaciers to both long- and short-term climate variability is
crucial for effective environmental resource planning and hazard mitigation (Bonsal et al., 2020;
Moore et al., 2009). Water resource management must adapt in communities downstream from
alpine glaciers to larger influxes of winter and early spring streamflow and account for long-term
implications of water availability (Marshall, 2011; Menounos et al., 2020).
1.1 Motivations and objectives
The primary objective of my research is to contribute to our understanding of past glacier
fluctuations in western Canada to help inform broader-scale climatic interpretations (Balco,
2020a; Menounos et al., 2009). To accomplish this, I developed a new glacier chronology for
glaciers among the Itsi Range in the Selwyn Mountains of the northeastern Canadian Cordillera.
Additionally, I aim to address the methodological uncertainties of geologic scatter that underpin
TCND research in this region. My research seeks to answer the following questions:
1. How have western Canadian glaciers responded to late Quaternary climate change?
2. Is there a plausible cause for observed geologic scatter in this region apart from TCN
inheritance or post-depositional disturbance?
3. How does the newly derived glacier chronology correlate with existing moraine
records, regional advances, and climate proxy data to enhance our understanding of
regional Holocene climate drivers?
1.2 Document structure
Chapter 2 presents new numerical ages of cirque glacier moraines in the northeastern
Canadian Cordillera using TCND methods and derive geodetic mass balance estimates among
glaciers of the Itsi Range since 1949. Chapter 3 compares the new numerical ages and positions

2

of moraines with the western Canada moraine record and discusses observations of geologic
scatter. Chapter 4 presents further analysis of the new moraine record with relevant climate
proxy records. Chapter 5 recounts the significance of my findings in the context of paleoclimate
research, explores remaining knowledge gaps, and the methodological advancements required to
improve the accuracy and correlative potential of TCND derived surface exposure ages.

2. Holocene activity among the Itsi Range
Chapter 2 considers how glaciers in the Itsi Range fluctuated during the late Holocene. I
developed my glacier chronology by using terrestrial cosmogenic 10Be surface exposure dating
techniques.
2.1 Climate and the cryosphere
2.1.1 Spatial and temporal climate scales
The cryosphere encompasses all frozen forms of water on Earth and is highly sensitive to
changes in mean global temperature (Marshall, 2011; Nye, 1960). Climate controls the
precipitation and temperature conditions that regulate glacier formation. Glaciers are perennial
masses of ice that flow due to gravitational deformation under their own weight (Cogley et al.,
2011; Marshall, 2011; Roe and O’Neal, 2009). Glaciers form at localities where snow
accumulation equals or exceeds ablation multi-annually, leading to the compression of snow into
dense glacier ice. Glacier retreat occurs when ablation exceeds accumulation (Marshall, 2011).
Mass balance refers to this equilibrium between accumulation and ablation (Cogley et al., 2011).
Holocene climate drivers operate at several different spatial and temporal scales within
western Canada; explaining the behaviour of glacier fluctuations requires considering these
scales (Clague et al., 2009; Menounos et al., 2009). Three dominant frequencies of climate
3

change characterize the Holocene: millennial, decadal-to-centennial, and sub-decadal (Clague et
al., 2009). The reduction of Northern Hemisphere (NH) summer solar insolation is responsible
for millennial-scale periods of successive glacier expansion throughout the Holocene (Clague et
al., 2009; Wanner et al., 2011) (Figure 2.1). Changes in solar activity, volcanism, and
teleconnection patterns drive decadal-to-centennial glacier fluctuations (Porter, 1986; Wiles et
al., 2004). Non-climatic factors such as topography, elevation, and glacier hypsometry affect
sub-decadal responses (Clague et al., 2009; Menounos et al., 2013b; Wanner et al., 2011).

Figure 2.1: July solar insolation at 65° N (W/m2) (Berger and Loutre, 1991) and NH 50 year
ensemble mean temperature anomalies (Erb et al., 2022) throughout the Holocene.
Additionally, climate broadly operates on three spatial scales: planetary, continental, and
regional. On a planetary scale, incident radiation from the sun reaches Earth's surface at lower
angles towards the poles, resulting in a cooler climate at higher latitudes from the reduced heat
flux and vice versa at lower latitudes (Berger, 1988). At a continental scale, ocean temperatures
strongly influence land temperatures due to water’s high capacity for storing solar energy
(Dommenget, 2009; Menounos et al., 2020). Topography affects climate at a local-to-regional
scale; the troposphere cools with increases of elevation, and rain shadows cause desert climates
on the leeward side of mountain ranges (Clague et al., 2009).
4

2.1.2 Glacier chronology as climate proxy
One of the greatest resources for understanding the natural variability of past climate is
reconstructing past glacier behavior through the development of glacier chronologies (Balco,
2020a, 2011; Menounos et al., 2009). Mountain glacier mass balance in the NH is most sensitive
to variations in boreal summer (June through August) temperatures and, to a lesser extent
precipitation, responding with geometric adjustments of thickness, length, and width (Balco,
2020a; Cogley et al., 2011; Hodge et al., 1998; Menounos et al., 2009; Nye, 1960; Solomina et
al., 2015). Due to lateral constraints, many valley and cirque glaciers exhibit fluctuations of
length as their most common responses to long-term changes in climate (Menounos et al., 2009).
Although long records of glacier mass balance offer the most direct way to infer how climate
affects glacier health, ample evidence of glacier length change exists and offers the ability to
glean information about climatic change (Oerlemans, 1994).
Non-climatic effects such as response time or how glaciers may respond to random, highfrequency weather events complicate insight into past climate (Roe and O’Neal, 2009). Glacier
response, in terms of fluctuating geometry, to changes in climate lags at a timescale ranging from
several years to centuries (Jóhannesson et al., 1989). This response may indicate either the
natural interannual variability of a stationary climate (i.e., one that contains a constant mean and
constant standard deviation), which does not reflect a true change in climate forcing, or a
statistically significant change in the expected response to natural climate variability (Roe and
O’Neal, 2009).
The null hypothesis, that any change in glacier length must be temporally and spatially
greater than those expected to be driven by the observed record of interannual variability in a
constant climate, must be rejected to demonstrate a glacier chronology as explicitly driven by
5

climate change (Roe and O’Neal, 2009). The signal of anthropogenic warming evidenced by a
globally synchronous retreat of glaciers since the late-20th century has been demonstrated as
statistically significant, exceeding the natural interannual variability of the climate system;
however, even the LIA advances cannot be demonstrated to exceed this natural variability
(Reichert et al., 2002; Roe and O’Neal, 2009). Roe and O'Neal (2009) advocate that lateHolocene moraine chronologies merely represent the random noise of a complex climate system
and that the inherently stochastic glacier response to climate prevents asserting any significant
correlations with paleoclimate drivers (Kirkbride and Winkler, 2012; Roe and O’Neal, 2009).
Although this perspective may be provocative and circumspect, the criticism offers a valid
rationale for examining quantified glacier response time lags before interpretation of dated
terminal and recessional moraines.
Bracketing ages of a given glacier advance remains an important prerequisite to closely
link climate drivers with glacier response as geologically evident glacier fluctuations can
preserve decadal, centennial, and millennial scale climatic signals indicating regional and, when
viewed synoptically, global Holocene climate trends (Hawkins et al., 2021; Kaplan et al., 2016;
Menounos et al., 2009; Schaefer et al., 2009; Schimmelpfennig et al., 2022; Solomina et al.,
2015).
Glacier response time operates as a function of mean thickness, terminus ablation rate,
hypsometry, and mass balance gradient, therefore highly variable among different environments
and geomorphological settings (e.g., high-relief coastal and low-relief continental interior
glaciers will have respectively fast and slow response times) (Kirkbride and Winkler, 2012). This
variability of response time implies that individual glaciers will respond differently, as
interpreted by their geomorphic record, to the same climatic history, rendering moraine
6

formation scale-dependent wherein large glaciers filter high-frequency climate signals (i.e.,
overriding of previously deposited moraines) and small glaciers respond more readily (Kirkbride
and Winkler, 2012).
2.1.3 Modern cryospheric conditions
Currently, perennial ice covers 10.8% of the Earth’s land surface with an additional
15.4% underlain by permafrost – frozen stores of soil at the Earth's veneer (Marshall, 2011). Net
radiation from the sun governs temperatures on Earth; however, because the change in solar
radiation on Earth over time is quite small, the greatest contributing factor to planetary climate is
global albedo (i.e., the planet's reflectivity) which largely depends upon snow and ice cover
(Marshall, 2011). Due to its high albedo, the cryosphere reflects shortwave radiation from the
sun at high latitudes and mountainous regions back into space which helps keep Earth's surface
temperature stable (Groisman et al., 1994). This albedo reduces incoming net solar radiation
which, with sufficient cryospheric cover, contributes to a reduced global mean annual
temperature (Groisman et al., 1994; Marshall, 2011). Cryospheric albedo also affects regional-tolocal climate due to its influence on atmospheric circulation (Marshall, 2011).
Global temperature variability on the order of 1 degree Celsius has driven Holocene
fluctuations of alpine glaciers (Marshall, 2011). Recent warming since the late 19th century is not
in line with expected orbital cooling trends (IPCC, 2021; Kaufman and Broadman, 2023).
Increased atmospheric Greenhouse Gas (GHG) emissions, as well as solar and volcanic activity,
have primarily driven trend divergence (IPCC, 2021; Solomina et al., 2015). Arctic warming
increases at over two times the average rate of global surface temperature warming, making this
a vulnerable region to climate warming (IPCC, 2021; Marshall, 2011).

7

With unprecedented increases in GHG emissions since the late 20th century, tropospheric
warming has destabilized the cryosphere, impacting hemispheric circulation, sea level rise, and
ice-albedo feedbacks leading to a warming climate (IPCC, 2021). Additionally, the potential for
a dramatic increase of carbon emissions from permafrost stores has risen as permafrost melt
leads to a sustained release of GHGs CO2 and CH4 due to positive feedback effects from climate
warming (IPCC, 2021; Marshall, 2011). Onset of anthropogenic climate change, resulting in the
global reduction of glacier mass, diminishes the capacity for glacier and snow melt to support
regional energy demands, ecosystems, fisheries, and agriculture of western Canada and increases
the likelihood of associated glacier hazards like outburst floods (IPCC, 2022; Milner et al., 2017;
Rounce et al., 2023). Global projections reveal that glaciers will experience mass loss ranging
from 28 ± 9% (SSP1-2.6) to 44 ± 20% (SSP5-8.5) by the year 2100 relative to 2015 (Rounce et
al., 2023).
2.1.4 Canadian Cordillera cryospheric activity
The most recent glaciation, the Last Glacial Maximum (LGM), occurred from 23 – 19 ka
and covered a significant portion of the NH with ice sheets (IPCC, 2021). The Cordilleran and
Laurentide Ice Sheets (CIS/LIS) merged around 22 ka, covering much of modern Canada and the
northern United States (Dalton et al., 2023; Darvill et al., 2018). The CIS attained its maximum
extent circa 18.5 ka with subsequent detachment and retreat of the two ice sheets beginning at
17.5 ka (Dalton et al., 2023). Following the LGM towards the end of the Pleistocene, global
climate fluctuated between glacial and interglacial conditions (Menounos et al., 2017). A release
of warm deep-ocean waters in the North Atlantic that had built up over time reignited the
Atlantic meridional overturning circulation causing the Bølling-Allerød (BA), an abrupt
warming event at 14.6 ka (Thiagarajan et al., 2014). BA warming induced CIS mass thinning
8

which subsequently allowed valley and cirque moraines in the western Canadian Cordillera to
advance into newly deglaciated alpine areas as early as 14.0 ka (Menounos et al., 2017). The BA
represents the final interstadial of the last glacial period which terminated at 12.9 ka with the
onset of the Younger Dryas (YD), the Pleistocene’s final stadial (Broecker et al., 2010). The CIS
disappeared by the end of the YD which terminated at 11.7 ka (Menounos et al., 2017;
Rasmussen et al., 2006). A return to consistent interglacial conditions after 11.7 ka marked the
beginning of the Holocene (Rasmussen et al., 2006). The Holocene Thermal Maximum (HTM)
was a period during the early Holocene from approximately 11 – 5 ka (Renssen et al., 2009;
Schimmelpfennig et al., 2022). The HTM signal is most strongly recorded in the NH,
characterized by a relatively warm climate due to orbitally forced high summer insolation
(Carlson et al., 2008; Renssen et al., 2009). Decreasing insolation in the NH during the boreal
summer led to Neoglaciation, characterized as the period of time when NH glaciers advanced
from early-to-mid Holocene, namely HTM, minimal extents and encompasses multiple smaller
scale phases of advance and retreat (Luckman et al., 1993; Menounos et al., 2009; Porter and
Denton, 1967; Wanner et al., 2011). This reduction of NH insolation culminated as globally
recorded negative temperature anomalies that characterize the most recent cold relapse, the LIA,
occurring between 0.7 – 0.15 ka (1200 – 1800 CE) (Berger, 1978; Wanner et al., 2011).
Glaciers within the Canadian Cordillera expanded throughout the Holocene with
maximum downvalley positions of glaciers achieved in the last millennium during the LIA;
because of this, scant morphological evidence for previous Neoglacial activity in the Canadian
Cordillera remains as prior glacier activity, in the form of lateral and terminal moraines, is
commonly overridden and erased through a process of erosion censoring called “obliterative
overlap” (Gibbons et al., 1984; Kirkbride and Winkler, 2012; Luckman et al., 1993; Menounos et
9

al., 2009; Reyes et al., 2006). Erosion censoring by means of obliterative overlap effectively
erases former records of climate among geomorphological deposits. However, some highelevation pre-LIA cirque and valley moraines in the Yukon Territory (YT) and Northwest
Territories (NWT) were discovered hundreds of meters to several kilometers beyond LIA
extents, either forming shortly after or interacting with the decaying CIS (Lakeman et al., 2008;
Menounos et al., 2017). Climatically, anthropogenic warming has defined the last two centuries,
largely due to the introduction of GHGs into the atmosphere and has prematurely terminated
influence of Little Ice Age (LIA) climate (IPCC, 2021; Wanner et al., 2011). Compared to
baseline variability of centennial to millennial glacier fluctuations recorded in Holocene glacier
chronologies, the anomalous retreat of glaciers throughout the last century in alpine regions
worldwide provides evidence for anthropogenic climate change (IPCC, 2013; Solomina et al.,
2015).
2.2 Terrestrial in situ cosmogenic nuclide dating theory
Most of the Galactic Cosmic Radiation (GCR) that contributes to TCN production
originates from the Milky Way galaxy (Gosse and Phillips, 2001). While the GCR cascades from
its primary source to the Earth, interactions between primary GCR and nuclei of atmospheric
atoms create a series of reactions and particles with net energy eventually lost to the materials of
Earth, such as rock (Darvill, 2013; Davis Jr and Schaeffer, 1955). Production of cosmogenic
nuclides occurs first as primary cosmic rays, then as secondary cosmic rays, and finally as
nuclides produced through nucleonic spallation and muogenic reactions on certain elements
(Darvill, 2013; Schaefer et al., 2022). Various spatial and temporal factors such as latitude,
altitude, GCR intensity and solar modulation, atmospheric shielding, and the shifting
morphology of the landform surface through time determine the production rate of cosmogenic
10

nuclides in rocks exposed at the surface generated by this cosmic radiation flux (Darvill, 2013;
Lal, 1991; Schaefer et al., 2022). Terrestrial in situ cosmogenic nuclide (TCN) production rates
assume a reference condition that consists of a horizontal surface with no obstructions at high
latitude and at sea level. In practice, TCN samples collected for research cover a wide range of
latitudes, elevations, and surface slopes that may be partially obstructed. To accurately record
TCN exposure ages, implementation of a scaling model allows the determination of spallogenic
and muonic nuclide production rates for any latitude, altitude, depth, and time period (Lal, 1991;
Lal et al., 1960; Lifton et al., 2014).
After the discovery of cosmic radiation by physicist Victor Hess in 1912, geologists
promptly recognized the potential utility for solving geological problems with the first attempt at
TCN Dating (TCND) carried out in 1955 (Davis Jr and Schaeffer, 1955; Gosse and Phillips,
2001). When Accelerator Mass Spectrometry (AMS) was introduced in the 1980s, measurement
of cosmogenic nuclides became broadly accessible, thus applications of TCND grew
exponentially in the years to follow (Brown et al., 1982; Gosse and Phillips, 2001).
TCND has distinct advantages over other methods used to construct a chronology of
glacier fluctuations prior to the 20th century. Other approaches such as lichenometry,
dendrochronology, and radiocarbon dating of subfossil wood, rely on biological proxies.
Lichenometry, often criticized as pseudoscience, has been scrutinized for producing unreliable
ages and for lacking consistent methodological practices (Osborn et al., 2015). Although
radiocarbon dating can be highly valuable, it requires environments that support vegetative
growth and can only bracket an age range for glacier advance or retreat (Darvill, 2013; Gosse
and Phillips, 2001). The geochemical method of tephrochronology requires the presence of

11

volcanic tephra deposits within moraine stratigraphy making it less suitable for late Quaternary
research (Osborn et al., 2001; Preece et al., 2011).
The relevance of TCND for constraining glacier chronologies is based upon the
assumption that fresh rock surfaces are created at the bed of glaciers which are concealed from
exposure to the cosmic-ray flux (Balco, 2020a, 2011). Boulders eroded at the glacier bed are
carried to the ice margin where they are exposed during ice retreat or deposited as a moraine
(Balco, 2020a, 2011). Once emplaced, the rock becomes exposed to the cosmic-ray flux and
build-up of previous cosmogenic nuclide concentrations are considered negligible; therefore, the
duration of exposure for a relevant rock can be used as a direct proxy for the timing of glacier ice
stabilization and retreat (Balco, 2020a, 2011). Accurately dating a glacial landform requires
resolving two questions: 1) how long a sample has been exposed to the cosmic-ray flux and; 2)
whether or not the sample's exposure duration can be linked to the presence of a glacier at that
time (Balco, 2020a).
I focus on applications of 10Be for dating in this thesis due to its well constrained
production rates (within the range of 3 – 4% uncertainty), abundance in the ubiquitous mineral
quartz, smallest amount of uncertainty in terms of coefficients of variation, and ability for use at
the centennial-scale throughout the Holocene (Balco, 2020a, 2011; Gosse and Phillips, 2001;
Schaefer et al., 2022). Technical advances in the 21st century now allow analysis of LIA era 10Be
ages producing low uncertainties and exposure ages consistent with historical records (Braumann
et al., 2020; Schaefer et al., 2022, 2009).
Successfully carrying out TCND methods requires knowledge of TCN production at
depth below the Earth's surface. Production by spallation, which occurs within the first 3 m of
the Earth’s surface, produces 96.4% of 10Be; below 3 m the neutron flux diminishes to <1% and
12

muonic reactions dominate the remaining 3.6% of 10Be production (Darvill, 2013; Heisinger et
al., 2002; Hippe, 2017). Quartz (SiO2) is the primary target for TCND samples in this research
context. High-energy neutron spallation on elements Si and O produce 10Be isotopes,
specifically: beryllium with two additional neutrons (10Be) from oxygen [16O(n,4p3n)10Be] and
from silicon [28Si(n,x)10Be] (Gosse and Phillips, 2001). Exposure time depends on the
concentration of 10Be:

𝑁10 =

𝑃10
[1 − 𝑒 −𝜆10 𝑡 ]
𝜆10

(1)

Where N10 is the 10Be concentration in quartz (atoms g-1), P10 is the production rate for 10Be
(atoms g-1 yr-1), 𝜆10 is the 10Be decay constant (4.99 x 10-7 yr-1), and t is the exposure time (years)
(Balco, 2011).
2.3 Study site
Evolution of the Canadian Cordillera began around 750 Ma with the opening of an
ancestral ocean basin due to rifting of the Neoproterozoic supercontinent Rodinia (Monger and
Price, 2002). A deep-water depositional environment formed at this time, known as the Selwyn
Basin tectonic assemblage, which accumulated clastic sediments from the Late Neoproterozoic
to the Triassic (Gordey, 2013; Smith and Bennett, 2004). Beginning with the Jurassic period,
western offshore subduction zones converged with the North American continent wherein
processes of accretion formed a new continental margin by the Late Cretaceous (Monger and
Price, 2002). During the Mid-Cretaceous, a continental arc formed the Selwyn Plutonic Suite that
intruded into the sedimentary deposits of the Selwyn Basin forming quartz-rich granitic massifs
that host the alpine glacier and genetically related Neoglacial moraine deposits sampled for this
research (Gordey, 2013; Smith and Bennett, 2004).
13

The Selwyn Mountains straddle the YT and NWT border from 61°N to 64°N and belong
to the Taiga Cordillera ecozone (Smith and Bennett, 2004). High elevation peaks separated by
broad valleys characterize the region. The highest point in this mountain range is Keele Peak at
2970 m, with local relief ranging between 900 – 1500 m. Mid-Cretaceous (98 – 93 Ma) quartz
monzonite of the Mayo Suite characterize the dominant bedrock geology of the Itsi Range,
matching the lithology of the moraine boulders and erratics chosen for sampling (Colpron,
2022).
Buried thrust faults encompass much of the Selwyn Mountains and the geographically
adjacent Mackenzie Mountains have been documented to experience active thrust and strike-slip
faulting generating high-magnitude earthquakes (Gordey, 2013; Hyndman et al., 2005). The
earthquakes recorded in these northeastern Canadian Cordillera mountain fronts trigger
rockslides along arêtes and cirque glacier headwalls (Smith and Bennett, 2004). Within a 100 km
radius of the Itsi Range, 165 seismic events have been recorded by the National Earthquake
DataBase (NEDB) since 1985; of these earthquakes, four range in magnitude from 4.6 to 5.2.
The Selwyn Mountains contain the highest values of precipitation outside of YT’s coastal ranges.
Among the Pelly Mountains just southwest of the Selwyn Mountains, rock avalanches generating
volumes of non-carbonate rock talus on the order of 1 to 5 x 106 m3 have been reported (Jackson
and Isobe, 1990).
Itsi Mountain (~2500 m a.s.l.) and its neighbouring landforms comprise the Itsi Range, a
massif located within the Selwyn Mountains of the eastern YT. The Itsi Range covers an area of
approximately 150 km2, corresponding to locally exposed contiguous pre-Quaternary bedrock
(Jackson et al., 1993). The Itsi Range contains evidence of Quaternary glaciation through the
presence of cirque glaciers, arêtes, Neoglacial moraines, and expansive U-shaped valleys.
14

Assessment of August 2022 Sentinel-2 imagery shows 26 active cirque and niche glaciers
comprising a summed total area of 11.30 km2 (Figure 2.2). Northern slope-facing cirque glaciers
have retained a greater volume of ice than their southern-slope facing counterparts. Mean Annual
air Temperature (MAT) and Mean Annual Precipitation (MAP) for the Itsi Range are
interpolated from a raster grid as -8.5 °C and 519 mm yr-1, respectively (Fick and Hijmans,
2017).
I completed an extensive review of potential study areas within the Selwyn and
Mackenzie Mountain ranges of the YT and NWT before selecting a site among the Itsi Range
(Figure 2.2). I first examined potential moraines using historic aerial stereoscopic imagery
available at the Yukon Energy, Mines, and Resources Library in Whitehorse, YT. This led to the
identification of a suitable pair of Holocene moraines and erratics down valley of the moraine
pair which sit on the forefield of an active unnamed south-facing cirque glacier, herein referred
to as the Itsi Glacier (Figure 2.4). The closely spaced moraines derive from quartz monzonite
bedrock and lie several hundred meters down valley from cirque headwalls. ArcticDEM reveals
older morphological evidence of the Itsi Glacier's activity, indicating the presence of latest
Pleistocene lateral valley moraines. The Itsi Glacier likely constructed these moraines as it
advanced soon after, or concurrently with, the decay of the CIS (Menounos et al., 2017).
The Itsi Glacier lies in an area of discontinuous permafrost (50 – 90% coverage; Figure
2.7) which places it in the middle of the basal thermal regime process-form continuum, with the
ice-margin likely sensitive to decadal or seasonal climate forcings (Evans, 2009; Obu et al.,
2019). The Itsi Massif lies in a setting of low- to moderate-relief, with an approximate vertical
distance between valley floor and summit ridge of 1250 m and an average relief of 163 m
between Itsi Glacier and its northern headwall.
15

16

Figure 2.2: (a) Regional view of the study site in relation to adjacent physiographic provinces.
(b) Itsi Range study area and glaciers overlayed on Sentinel 2 10 m true colour imagery captured
on 11 August 2022.
2.4 Methodology
2.4.1

10Be dating

Given the expense required to access the Itsi Range, funding for TCN analysis was
limited to 10 samples; an insufficient amount to attain a synoptic understanding of site-specific
paleoclimate dynamics (Balco, 2020a). Consequently, I ascertained that augmenting ongoing
research in the eastern YT/western NWT to establish a more robust regional chronology would
be more beneficial than simply increasing sample numbers at previously studied sites.
We accessed the field site via a chartered Bell 407 helicopter during late summer (13
August 2022) to minimize the presence of seasonal snow cover. Based on criteria outlined by
Darvill (2013) and Gosse and Phillips (2001), we chose ten boulder samples in the field: five
boulders from the outer moraine, three boulders from the inner moraine, and two glacial erratics.
Relevant qualitative site selection criteria include identifying a quartz-rich lithology, selecting
the tallest boulders relative to the morainal ridge and closest to the moraine crest (avoids postdepositional shielding), selecting large boulders indicative of stable emplacement within the
moraine matrix (avoids post-depositional exhumation), and ensuring sampled surfaces are
sufficiently flat and horizontal showing minimal signs of exhumation and erosion (e.g., presence
of glacial striations and a moderate amount of lichen cover present).
Tomkins et al. (2021) make a case that no conclusive pattern to the presence of outliers in
a glacier chronology exists when comparing those sampled from three intra-landform locations:
moraine crests, ice-proximal slopes, and ice-distal slopes. The most important geomorphic
consideration to mitigate geologic scatter remains moraine sedimentology, with boulder-rich
17

matrix-poor moraines more likely to yield reliably accurate boulder ages than matrix-rich
unconsolidated moraines (Tomkins et al., 2021). We followed recommendations from prior
literature and practice by sampling from moraine crests as Itsi Glacier moraines exhibited a
boulder-rich and matrix-poor composition (Gosse and Phillips, 2001; Hawkins et al., 2021).
We used a gas-powered rock saw to recover samples for TCN analysis. We made long,
straight incisions along the length and width of the chosen surface, creating a grid of extractable
rock ‘brownies’ with dimensions roughly equal to 5 x 5 x 3 cm (L x W x H) pried out using a
hammer and chisel. The process resulted in approximately ~1 kg of extractable material per
boulder with samples stored in cloth bags labeled with site identifier, year, and sample number
(Figure A.4). I recorded the following properties for each boulder sample: sample ID, strike/dip,
boulder dimensions, lithology, amount of vegetation cover, and geomorphic context (Table 2.1).
At each boulder’s location, a Magellan eXplorist 210 GPS unit recorded latitude, longitude, and
elevation. A Brunton transit compass with clinometer recorded the topography of the landscape
at 45° intervals to account for topographic shielding of cosmic rays for each boulder. A roughly
15 x 20 cm surface of freshly cut rock remains after extracting each sample. We placed small
rocks atop the modified boulder surface to minimize aesthetic disturbance which risks
compromising the site’s cultural and spiritual values (Schaefer et al., 2022).
Table 2.1: Itsi Glacier field data collected on 13 August 2022.
Sample
Latitude
ID
(d. d.)
Outer moraine
228B01
62.8904
228B02
62.8900
228B03
62.8901
228B04
62.8927
228B05
62.8928
Inner moraine

Longitude
(d. d.)

Elevation
(m)

Strike°/
Dip°a

Shielding
Factorb

L
(m)

W
(m)

H
(m)

-130.2817
-130.2803
-130.2799
-130.2905
-130.2903

1722
1724
1722
1750
1761

108/02
342/08
214/04
138/07
104/04

0.983
0.991
0.991
0.986
0.986

1.45
1.80
3.18
1.50
1.25

0.85
1.60
2.24
1.10
0.75

0.70
0.80
1.60
0.55
0.60

18

228B06
228B07
228B08
Erratics
228B09
228B10

62.8938
62.8938
62.8936

-130.2913
-130.2910
-130.2916

1796
1798
1789

160/12
290/15
260/13

0.987
0.986
0.986

2.40
2.20
2.40

2.20
1.00
2.10

0.90
0.70
1.10

62.8921
62.8921

-130.2916
-130.2916

1757
1757

014/05
326/59

0.982
0.821

1.50
2.05

1.45
1.00

0.65
1.00

a

Strike and dip of sampled surface (right-hand rule); values recorded with magnetic declination set to 0 degrees.
World Magnetic Model (WMM) magnetic declination for the site is 19.73° E ± 0.55°.
b
Shielding factor calculated using the topographic shielding calculator found at https://stoneage.iced.org/math/skyline/skyline_in.html (version info.: wrapper: 2.0; skyline: 2.0).

Sample 228B10, an erratic boulder immediately adjacent to 228B09 and displaying the
same highly weathered and lichen covered characteristics, had its geometry sampled
experimentally to see how an erratic displaying poor geomorphic characteristics in the context of
ideal 10Be sampling would be affected when compared to 228B09, providing insight into the
suitability of standard field sampling recommendations. 228B10’s surface strike/dip of 326/59
and shielding factor of 0.821 augments the sample’s geometric complexity, foreshortening its
exposure to the cosmic-ray flux, i.e., abating the uniform cosmic-ray exposure history that a
sufficiently flat and horizontal surface would have (Gosse and Phillips, 2001).
Upon collection, samples were sent to the Lamont-Doherty Earth Observatory (LDEO) in
Palisades, New York for quartz extraction. Before the samples are subjected to the physical
processes of crushing and sieving, they must be described in terms of average thickness and
weight. A weighted average thickness of all pieces of an individual sample utilized for extracting
quartz is then calculated and recorded. Samples were crushed using a hydraulic press, milled,
and sieved to fractions varying between 63 and 710 μm. The rock crushing lab was thoroughly
cleaned between samples to avoid 10Be contamination. We attempted the LDEO’s (2012a)
standard methods of quartz extraction from whole rock including techniques of froth flotation,
heavy liquid mineral separation, and up to five leaches of 5% hydrofluoric acid (HF) per sample,
but this was unsuccessful as the rock samples completely dissolved in HF yielding no extractable
19

quartz (LDEO, 2012a). Subsequent XRD analysis of samples 228B01, 228B07, and 228B09 at
the UNBC Northern Analytical Laboratory Services revealed a substantial percentage of quartz
present in each rock sample; consequently, Rinterknecht’s (2003) alternative procedures of Ophosphoric boiling were followed successfully. Sample 228B02 was not processed due to a lack
of available material and time constraints.
Young samples of rock for use in cosmogenic surface exposure dating require
significantly more quartz for 10Be measurement than older samples. As our sampled late
Holocene moraines are particularly young, it was recommended to carry out an experimental
procedure of sample preparation. All samples were spiked with an average of 107.13 μg 9Be and
the remaining sample bulked with Fe instead of the typical ~250 μg pure 9Be carrier. Using less
9

Be carrier improves 10Be data precision when determining younger exposure ages. This

approach allows for improved accuracy of ages, counting statistics, and requires ~40% less
quartz compared to current standards (Braumann et al., 2021; Walker et al., 2022). Two low
concentration process blanks were included for the nine samples (Hippe, 2017; LDEO, 2012b).
Standard methods of Be extraction from quartz were used after 9Be spiking. This includes anion
exchange columns for removal of FeIII and cation exchange columns for the separation of Be
from Al and Ti (LDEO, 2012b). Sample aliquots of extracted Be were shipped to Lawrence
Livermore National Laboratory Center for Accelerator Mass Spectrometry (LLNL-CAMS) to
carry out AMS measurement for beryllium-isotope ratio analyses relative to the 10Be/9Be
standard 07KNSTD dilution series (Nishiizumi et al., 2007).
Surface exposure ages are reported and discussed as either years or thousands of years
(ka) Before Present (BP) using 1950 CE as the reference year. Exposure ages were calculated
using the online exposure age calculator formerly known as CRONUS-Earth v3
20

(https://hess.ess.washington.edu/) (Phillips et al., 2016). Ages are calculated using the LiftonSato-Dunai nuclide-dependent (LSDn) scaling framework with default Be atom reference
production rates as no regional 10Be production rate calibration site for the YT is available
(Borchers et al., 2016; Lifton, 2016; Lifton et al., 2014). Moraine ages and erratics are
conservatively reported as the median ± the interquartile range (IQR) rounded to the nearest
decade. The median ± IQR is less sensitive to extreme values and accounts for uncertainties such
as geological scatter, providing a more conservative estimate of the moraine's true age (Darvill et
al., 2022; Menounos et al., 2017). Measurement uncertainties for exposure age uncertainty
(Table 2.4) is calculated as propagated error:
√(𝜎𝑖𝑛𝑡 )2 + (𝜎𝑒𝑥𝑡 )2

(2)

Where 𝜎𝑖𝑛𝑡 is internal error and 𝜎𝑒𝑥𝑡 is external error.
2.4.2 DEM creation and geodetic mass balance
Unlike geometric fluctuations in glacier length, which reflect long-term changes, mass
balance records provide a more immediate response to changes in atmospheric conditions
(Oerlemans, 1994; Zemp et al., 2009). Geodetic mass balance studies reflect changes in surface
elevation based on assumed density values for water and ice (Beedle et al., 2014). To gain
insight into glacier response to climate change over the past century, I estimated the geodetic
mass balance for the Itsi Range using photogrammetric methods. To achieve this, I ordered aerial
photo negative scans and associated fiducial calibration reports from the Natural Resources
Canada (NRC) National Air Photo Library (NAPL) scanned at 2032 dpi and saved as TIF files
(Table 2.2). These air photos provide nearly full spatial coverage of the Itsi Range on 8/24/1949
and 9/9/2004.

21

Table 2.2: Aerial photographs used for glacier change analysis.
Year

Date

Roll
ID

1949

24 Aug.

A12255

2004

09 Sep.

A28519

Photo
Numbers
0367 - 0370,
0433 - 0436
0145 - 0148

Scale

GCPsa

RMSEb (m)
XY

RMSE (m)
Z

1:40 000

34

4.780

5.853

1:40 000

20

3.640

4.156

a

Total number of Ground Control Points (GCPs) used.
b
Root Mean Square Error (RMSE) of GCP positions.

I created a ground control point geopackage in QGIS (v. 3.36.1) to georectify the aerial
photos. I placed GCPs on stable, non-glaciated terrain using both a spaceborne DEM
(ArcticDEM) and Planet imagery, both resampled to a 4 m resolution. The ArcticDEM mosaic is
comprised of DEM strips taken between 9/12/2011 and 2/12/2021 (Porter et al., 2023). I aimed
for 5 – 10 GCPs per photo, placing them on consistent morphological features among the suite of
images such as pyramidal peaks, moraines, and distinctive bedrock features (Table 2.2). I saved
GCP coordinates to CRS EPSG:4326 (WGS84) and exported them as a CSV file.
I first optimized the scanned air photos in Python using the Contrast Limited Adaptive
Histogram Equalization (CLAHE) method from the OpenCV package to correct for film
overexposure on glacier ice (Reza, 2004). I used Agisoft Metashape Professional Edition (v.
2.1.2) for DEM and orthomosaic construction from the air photos. I processed each roll of air
photos separately, beginning with automated masking of scanned film borders and fiducial mark
detection. Fiducial coordinates from the calibration report were entered in the camera calibration
dialog box, with transformation type set to conformal. In the reference panel, the CRS was set to
EPSG:4326 and marker accuracy set to 20 m. To generate a sparse cloud of tie points, I then
aligned the photos at the lowest accuracy setting, set reference preselection to sequential, and
applied masks to key points. Lower accuracy settings coarsen resolution but minimize voids in
the DEM. After GCP import and assignment to each image, I built a point cloud at the lowest
accuracy setting. I created the DEM from this point cloud and projected it to EPSG:3579
22

(NAD83(CSRS)/Yukon Albers). I then constructed an orthomosaic of the aerial photos based on
the DEM (Figure 2.3).

Figure 2.3: Three-dimensional visualization of the 2004 orthomosaic of the Itsi Range, optimized
using CLAHE. No vertical exaggeration is applied. Glaciers are outlined in red, outer moraine
highlighted in yellow, and inner moraine highlighted in blue. Scale is approximate. Created in
QGIS using the Qgis2threejs plugin.
To assess geodetic mass balance over the periods 1949 – 2004 and 2004 – 2022, I
imported the newly georectified DEMs, the ArcticDEM mosaic, and an additional ArcticDEM
strip captured on 9/19/2022 into a Python environment (v. 3.11.9) with the package xDEM (v.
0.0.19) (contributors, 2024; Porter et al., 2022). I identified fourteen glaciers in the Itsi Range
present across all three rasters, and glacier outlines for each date were saved as vector polygon
geopackages. I reprojected vectors and rasters to EPSG:3579 and resampled rasters to the lowest
resolution DEM (photo A12255) to 8.11 m. I utilized the ArcticDEM mosaic as the reference for
coregistration. The coregistration process required two non glacierized, stable ground inlier
23

masks created from 1949 and 2004 glacier extents. The aerial photo-derived DEMs then undergo
deramping by correcting directional biases using a 2nd order polynomial (Girod et al., 2017). I
further processed the deramped DEMs and the 2022 DEM using a nonlinear blockwise
coregistration approach, which divides each DEM into a 16 chunk grid (Nuth and Kääb, 2011). I
assessed normalized median absolute deviation (NMAD) errors before and after coregistration to
ensure an improved fit for each DEM (Table 2.3).
Table 2.3: NMAD errors before and after coregistration for the DEM of each year.
DEM year
1949
2004
2022

NMAD before (m)
10.61
4.70
3.65

NMAD after (m)
4.70
2.08
1.75

I created a difference-DEM (dDEM) for each period by subtracting the oldest year from
the newest, resulting in dDEM 2004 – 1949 and dDEM 2022 – 2004 (Figure A.6). Outlier
removal entailed using a distance filter for each dDEM with a radius of 1 pixel and mean outlier
threshold of 50. An additional filter removed elevation difference values ≥ 200 for dDEM 2004 –
1949 and values ≥ 50 for dDEM 2022 – 2004. I divided the filtered dDEMs into 50 m mean
elevation bins and fit to a third-order polynomial (Mukherjee et al., 2023). I used these fitted bins
to derive geodetic mass balance (𝐵) values in metres of water equivalent (m w.e.) using equation
(3:
𝐵=

𝑑𝑉 𝜌𝑖
⋅
𝐴 𝜌𝑤

(3)

Where 𝑑𝑉 is total volume change, 𝐴 is average area of glacier extents, and an assumed ice
density conversion factor (𝜌𝑖 ) of 850 kg m-3 and water density (𝜌𝑤 ) of 1000 kg m-3 (Huss, 2013;
Mukherjee et al., 2023).

24

Geodetic mass balance uncertainty is standardized using stable terrain as error proxy
(Hugonnet et al., 2022). I calculated combined errors (𝐸Δ ), such as area and elevation change
errors, by combining error terms (𝐸n ) as shown in equation (4 (Tennant and Menounos, 2013).
𝐸Δ = √𝐸12 + 𝐸22 +⋅⋅⋅ +𝐸𝑛2

(4)

Following methods by Menounos et al. (2019), I calculated random uncertainty on
geodetic mass balance (σΔ𝑀,𝑟𝑑𝑛 ) which derives from three independent sources of error: rate of
elevation change uncertainty (𝜎Δ𝑧 ), glacierized area uncertainty (𝜎𝐴 ), and volume to mass
conversion uncertainty (𝜎𝑓Δ𝑉 ). After correcting for effective sample size, I derived elevation
change uncertainty (𝜎Δ𝑧 ) from uncertainty on stable terrain (σΔℎ ):
σΔ𝑧 = σΔℎ ⋅ √

𝐴𝑐𝑜𝑟
5𝐴

(5)

Where 𝐴 is glacier area and 𝐴𝑐𝑜𝑟 = 𝜋𝐿2 . I assume a decorrelation length (𝐿) of 500 m
(Menounos et al., 2019). I derive glacierized area uncertainty (σ𝐴 ) using the summed perimeter
of each glacier (𝑃) and half the resolution of a pixel (𝑑𝑥). I set 𝑑𝑥 to match the resampled raster
resolution of 8.11 m.
𝑑𝑥
σ𝐴 = 𝑃 ⋅ ( )
2

(6)

I use these uncertainties to calculate non-systematic volume change uncertainty (σΔ𝑉 ):
2

σΔ𝑉 = √(σΔ𝑧 ⋅ (𝑝𝐴 + 5(1 − 𝑝𝐴 ))𝐴) + (σ𝐴 ⋅ Δ𝑧)2

(7)

25

Where 𝑝𝐴 represents the percentage of glacierized area after accounting for voids in the dDEM
due to DEM construction and outlier filtering, and 𝐴 is glacier area. Random uncertainty on
geodetic mass balance (σΔ𝑀,𝑟𝑑𝑛 ) is then derived using volume change uncertainty (σΔ𝑉 ):
σΔ𝑀,𝑟𝑑𝑛 = √(σΔ𝑉 ⋅ 𝑓Δ𝑉 )2 + (σ𝑓Δ𝑉 ⋅ Δ𝑉)

2

(8)

Where Δ𝑉 is volume change, and the mass to volume conversion factor (𝑓Δ𝑉 ) and associated
uncertainty (σ𝑓Δ𝑉 ) are 850 kg m-3 and 60 kg m-3, respectively.
I calculated volumetric change error (𝐸Δ𝑉 ) following methods by Tennant and Menounos
(2013):
2
𝐸Δ𝑉 = √Δ𝐻 2 𝐸𝐴2 + 𝐴2 𝐸Δ𝐻

(9)

Where Δ𝐻 is mean elevation change, 𝐸ΔH is elevation change error, 𝐴 is glacier area used for
volume change calculations, and 𝐸𝐴 is glacier area error.
2.4.3 Termini tracking
I used the Glacier Termini Change Tracking (GTT) toolbox, developed by Urbanski
(2018) for ArcGIS, to calculate retreat rates of the Itsi Range glaciers using historic aerial
photographs and satellite imagery. Following aforementioned methods, I constructed a highresolution DEM (1 m) and orthomosaic from 9/9/2004 aerial photos which facilitated the precise
delineation of moraines among the Itsi Range. I traced successive glacial termini for each glacier
and date, storing them as 20 separate vector line shapefiles for estimating rates of glacial retreat.
The Python-based GTT toolbox connected lines between each glacier terminus and moraine
terminus proxy at 20 m intervals (Figure A.3). By using suitable nested moraines as proxies for

26

glacier termini and incorporating median moraine ages from the Itsi Glacier 10Be chronology, I
extended the timeline for analyzing termini retreat rates to the LIA.

27

2.5 Results
2.5.1 Itsi Glacier 10Be chronology

28

Figure 2.4: Itsi Glacier and moraines. Aerial photo and Sentinel-2 derived glacier extents are
overlayed for 1949, 2004, and 2022. DEM basemap was generated using stereoscopic aerial
image pairs from 2004. Individually 10Be dated boulder ages from each moraine are included and
reported with propagated error as LSDn exposure age BP.
A common geomorphologic feature crops up throughout the Mackenzie and Selwyn
mountains: the preservation of nested moraines in close proximity to one another among alpine
cirque glaciers (Figure 2.8) (Hawkins et al., 2023). The Itsi Glacier forefield contains two such
nested moraines: an outer and inner moraine with median surface exposure ages of 960 years ±
360 (n = 4) and 200 years ± 50 (n = 3), respectively (Figure 2.5; Table 2.4). Both moraines are
unvegetated, have minimal lichen cover, and are distinctly sharp-crested ridges suggesting a late
Holocene age. Sampled inner moraine, outer moraine, and erratic boulders lie at average
elevations of 1782 m, 1776 m, and 1757 m, respectively. In total, the Itsi Range contains nine
pairs of morphologically similar nested moraines. The average distance to the headwall for all
Itsi Range inner and outer moraines is 1268 and 1338 m, with corresponding average elevations
of 1625 m and 1611 m.

29

Figure 2.5: 10Be ages from two neoglacial moraines emplaced by the Itsi Glacier reported in ka
BP. Grey curves show ages of individual boulders rendered as Gaussian distributions. The
coloured lines show the sum of each individual Gaussian distribution for the 10Be ages sorted by
morphostratigraphic position. Box-and-whisker plots corresponding to each moraine denote the
median ± IQR, minimum, and maximum surface exposure ages.
Erratics beyond the outer moraine date to a median exposure age of 11.61 ka ± 0.54
(Figure 2.6). Erratic 228B09, sampled utilizing standard methods, returned an exposure age of
10539 years ± 706. Erratic 228B10, sampled following the experimental method described
earlier, returned an exposure age of 12686 years ± 917. This yielded results that were the inverse
of what I expected as sampling 228B10 at a steep angle (59°) should foreshorten its exposure to
the cosmic-ray flux, reduce 10Be production by spallation, and consequently return a younger
age.

30

Figure 2.6: 10Be ages from two erratic boulders outboard the Itsi Glacier moraines reported in ka
BP. The pink line shows the sum of each individual Gaussian distribution. Box-and-whisker plot
denotes the median ± IQR, minimum, and maximum surface exposure ages of the two erratics.

31

Table 2.4: Itsi Glacier 10Be sample information.
Sample
ID

Latitude
(d.d.)

Longitude
(d.d.)

Elevation
(m a.s.l.)

Thickness

Shielding

(cm)

Quartz
(g)

Blankcorrected

Carrier
added

10

Be/9Be ratio

(g)a
Outer moraine
228B01
62.8904
228B03
62.8901
228B04
62.8927
228B05
62.8928

Inner moraine
228B06
62.8938
228B07
62.8938
228B08
62.8936
Erratics
228B09
228B10

62.8921
62.8921

-130.2817
-130.2799
-130.2905
-130.2903

-130.2913
-130.2910
-130.2916

-130.2916
-130.2916

1722
1722
1750
1761

1796
1798
1789

1757
1757

1.51
1.88
1.71
2.23

2.33
2.43
2.11

1.39
4.74

0.983
0.991
0.986
0.986

0.987
0.986
0.986

0.982
0.821

41.140
40.141
23.153
31.388

40.986
28.884
30.738

13.974
20.098

0.103
0.102
0.102
0.103

0.103
0.102
0.102

0.102
0.102

1.52E-13
1.04E-13
1.53E-13
4.98E-14

3.45E-14
1.68E-14
3.85E-14

4.42E-13
6.20E-13

Blankcorrected
10

Be/9Be ratio

Blankcorrected
10

Be conc.

Blank-corrected
10

Exposure

Exposure

AMS

Be conc. 1𝜎

age (LSDn)

age
uncertainty

facility

(yrs)b, c

(yrs)

1156
761
2094
432
Median ±
IQR

82
59
145
38

1𝜎 uncertainty

(atoms g−1)

uncertainty (atoms
g−1)

2.91E-15
2.63E-15
3.02E-15
1.55E-15

2.62E+04
1.82E+04
4.66E+04
1.10E+04

5.75E+02
5.08E+02
1.05E+03
3.82E+02

1.48E-15
1.35E-15
1.48E-15

8.42E-15
1.67E-14

5.72E+03
3.71E+03
8.55E+03

2.24E+05
2.18E+05

2.73E+02
3.46E+02
3.65E+02

4.84E+03
6.28E+03

196
109
317
Median ±
IQR
10539
12686
Median ±
IQR

LLNL-CAMS
LLNL-CAMS
LLNL-CAMS
LLNL-CAMS

960 ± 360
24
26
33

LLNL-CAMS
LLNL-CAMS
LLNL-CAMS

200 ± 50
706
917
11610 ±
540

LLNL-CAMS
LLNL-CAMS

a

Be carrier concentration for all samples was 1045.9 ug g-1.
Ages are reported as years BP (1950 CE) and were calculated using the online exposure age calculator formerly known as CRONUS-Earth v3 found at
https://hess.ess.washington.edu/ (last access: 13 November 2023) (wrapper 3.0.2, muons: 1A, constants as of 26 August 2020). Ages are calculated using the
Lifton-Sato-Dunai (LSDn) scaling framework with the default Beryllium atom production rate.
c
Moraine ages and erratics are reported as the median ± IQR rounded to the nearest decade.
b

32

2.5.2 Termini tracking
Results from the GTT tool record net retreat for all time periods (Table 2.5 and A.2).
Yearly retreat rates ranged from 0.07 ± 0.14 to 8.31 ± 11.15 m yr⁻¹. Notably, the period from
2004 to 2022 records a marked increase in yearly retreat rate compared to the earlier periods with
a more than two-fold increase from 1949 to 2022. The Itsi Glacier lost 50% of its glaciated area
between 1949 and 2022 and contributed to a 43% percent reduction in total area for Itsi Range
glaciers.
Table 2.5: Statistics derived from the GTT tool and area analysis among all glaciers of the Itsi
Range.

a

762

Mean yearly
retreat rate
(m yr-1)
0.07

Standard
deviation (m
yr-1)
0.14

1754 – 1949p

196

1.95

1949 – 2004o

55

2004 – 2022o

18

Time period
(CE)a

Elapsed
yearsb

992 – 1754p

p

Mean retreat
distance (m)

Glacier
countc

ΔArea
%

50

9

-

2.20

381

9

-

3.35

3.54

185

20

-26%

8.31

11.15

149

20

-23%

o

Time periods marked by and denote proxy and observed, respectively.
Elapsed years corresponding to time period.
c
9 of the 20 glaciers contained morphologically similar nested moraines appropriate for use as termini proxy. Two
glaciers have distinguished tongues with independently developed moraines and are assessed as separate glaciers for
retreat rate.
b

2.5.3 Geodetic mass balance
From 1949 to 2022, the total volumetric change among the 14 glaciers, which covered an
area of 16.307 km² in 1949, amounted to a loss of -0.269 km3 ± 0.207 (Table 2.6). The data
reveal a significant increase in the rate of mass loss in the latter period. From 1949 to 2004, the
glaciers record a negative geodetic mass balance of -0.190 m w.e. yr-1 ± 0.065. The period from
2004 to 2022, indicates the mass loss rate increased by 101%, with glaciers losing -0.382 m w.e.
yr-1 ± 0.136, denoting accelerated mass loss over the past two decades. Area change among all
glaciers of the Itsi Range follows similar trends (Table A.2).
33

Table 2.6: Mass budget of Itsi Range glaciers (1949 – 2022).
Start year
(CE)
1949
2004

End year
(CE)
2004
2022

Elapsed
Years
55
18

Analysis area
(km2)
14.448
11.180

ΔV
(km3)
-0.178 ± 0.173
-0.091 ± 0.034

Mass change
(m w.e. yr-1)
-0.190 ± 0.065
-0.382 ± 0.136

2.6 Discussion
2.6.1 Geologic scatter
Constraining and correlating millennial-scale Holocene climate events from exposuredated alpine glacier chronologies requires knowledge of geologic scatter for effective mitigation
(Balco, 2020a, 2011). Geologic scatter is a methodological issue stemming from geologic
processes concerning the scatter in exposure ages derived from glacial landforms (Balco, 2011;
Darvill, 2013; Heyman et al., 2011). Measurement uncertainty, post-depositional disturbance
(unconformably young ages), and inheritance (unconformably old ages) can cause excess
geologic scatter (Balco, 2020a, 2011; Heyman et al., 2011). Two primary issues are associated
with the occurrence of geologic scatter: 1) it affects the ability to estimate central tendency,
decreasing the accuracy of estimated sample ages and; 2) as geologic scatter reflects sample
dispersion, i.e., variance for a normal distribution, a high dispersion of ages from an individual
landform makes it difficult to attribute sample ages to a past climate event.
Averaging geological scatter produced by measurement uncertainty that is independent
with normal distributions will yield a more accurate age; however, taking the mean will not
improve accuracy of a calculated exposure age from scatter that results from geomorphic
processes which affects exposure of a sample to the cosmic rays necessary for TCN production
(Balco, 2011). Discerning between these two scatter types poses a challenge for TCND studies,
yet accounting for uncertainty related to geologic scatter permits the confident attribution of
34

derived ages as a direct result of climatic processes (Balco, 2020a, 2011; Goehring et al., 2022).
The degree of geologic scatter has been estimated for each moraine using the reduced chisquared (X2) statistic as defined by Balco (2011):
𝑋𝑅2 =

2
𝑛
1
𝑡̅𝑖
∑ [𝑡 −
]
𝑛 − 1 𝑖=1 𝑖 𝜎𝑡𝑖

(10)

Where 𝑛 is number of samples, 𝑡𝑖 is the apparent exposure age, 𝑡̅𝑖 is the mean exposure age, and
𝜎 is the measurement uncertainty (Crump et al., 2017). A calculated reduced X2 statistic of less
than 2 indicates minimal geologic scatter (Schaefer et al., 2009). Both moraines and erratics fail
to pass the chi-squared statistical test (Table 2.7).
Utilizing the X2 statistic assumes the sample distribution is normal (parametric), but this
may not be the case for exposure ages from moraine boulders due to geologic scatter. The
Kolmogorov-Smirnov (KS) test is a non-parametric test that can be used to determine if data are
significantly similar to or different from a normal distribution (Massey Jr, 1951). To carry out
the test, exposure ages are first standardized:
𝑡𝑖 − 𝑡̅
𝜎𝑡

(11)

Where 𝑡̅ is the mean and 𝜎𝑡 is the standard deviation. These standardized ages are then used for
the KS test:
𝐷𝑠𝑡𝑎𝑡 = 𝑚𝑎𝑥|𝐹𝑒𝑚𝑝𝑖𝑟𝑖𝑐𝑎𝑙 (𝑡𝑖 ) − 𝐹𝑛𝑜𝑟𝑚𝑎𝑙 (𝑡𝑖 )|

(12)

Where 𝐹𝑒𝑚𝑝𝑖𝑟𝑖𝑐𝑎𝑙 (𝑡𝑖 ) is the empirical cumulative distribution function (CDF) of the exposure
ages, 𝐹𝑛𝑜𝑟𝑚𝑎𝑙 (𝑡𝑖 ) is the CDF of the normal distribution, and 𝐷𝑠𝑡𝑎𝑡 is the maximum difference
between the two distributions (Massey Jr, 1951). The KS test is calculated in Python using the
35

SciPy package. However, the KS test may become unreliable when dealing with the small
sample sizes typical of moraine datasets (e.g., n = 2 to 6). Small sample sizes can return large pvalues that misrepresent conclusions about the data's distribution. Because of this, the KS test is
not used on individual moraine landform ages.
The non-parametric test Median Absolute Deviation (MAD) does not rely on the
assumption of normality, which may be a better fit when considering the issue of geologic scatter
skewing sampled moraine boulder distributions in this region. Unlike the mean and standard
deviation, which can be disproportionately influenced by outliers, MAD is calculated as the
median of the absolute deviations from the dataset's median, providing a better central tendency
estimate when dealing with non-parametric data (Rousseeuw and Croux, 1993).
(13)

𝑀𝐴𝐷 = 𝑚𝑒𝑑𝑖𝑎𝑛(|𝑡𝑖 − 𝑡̃|)

Where 𝑡𝑖 is the 𝑖 𝑡ℎ exposure age in the dataset and 𝑡̃ is the median exposure age of the dataset,
modified from Leys et al. (2013). A lower MAD value indicates that the data points are closely
clustered around the median, suggesting minimal scatter, while a higher MAD value points to
greater variability. In the context of exposure ages, a low MAD suggests that the ages are
consistent, indicating minimal influence from geomorphic processes, while a high MAD may
signal significant geologic scatter.
Table 2.7: Itsi Glacier moraine and erratic ages with accompanying reduced chi-squared (X2)
and Mean Absolute Deviation (MAD) statistics.
Position
Outer moraine
Inner moraine
Erratics

n
4
3
2

Median ± IQR
0.96 ka ± 0.36
0.20 ka ± 0.05
11.61 ka ± 0.54

X2
133.44
12.73
3.68

MAD
362
87
1073.5

Ice-cored moraines are moraines containing ice overlaid by sediment and disconnected
from the active glacial ice-margin (Lukas, 2011). An ice-cored moraine must be considered a
36

transitional landform as the establishment of a stable moraine requires the melt out of all ice
incorporated into the till (Lukas, 2011). The delay of surface stabilization caused by the melting
of dead-ice may not accurately reflect the age of the landform (Kirkbride and Winkler, 2012).
The thaw of ice-cored moraines has potentially contributed to observations of scatter (Goehring
et al., 2022; Hawkins et al., 2023) in this region as evidence for anthropogenic warming driving
geomorphic transformation of glacial deposits in the Canadian Arctic has emerged in the last
decade (Crump et al., 2017; Kokelj et al., 2017; Segal et al., 2016). As warming accelerates in
Arctic regions, the melting of permafrost in ice-cored moraines threatens to mobilize sediment
stores causing increased thaw slump activity (Kokelj et al., 2017; Obu et al., 2019; Segal et al.,
2016). Melt of dead-ice in ice-cored moraines, causing a thaw slump beneath the layer of
moraine debris, may contribute to post-depositional disturbance causing the dated exposure ages
to appear younger than the true landform age – a form of self-censoring (e.g., Crump et al., 2017;
Kirkbride and Winkler, 2012).

37

Figure 2.7: A regional overview of permafrost zone coverage in relation to the Itsi Range field
site from Obu et al. (2019).
The surficial geology of the field site is described as Neoglacial till wherein ‘end
moraines may be tens of metres thick and contain masses of buried glacial ice,’ lending credence
to the plausibility of thaw slump induced post-depositional disturbance (Jackson et al., 1993).
However, field observations of the moraines revealed that they are matrix-supported with little
evidence of permafrost thaw, such as thermokarst features or morphological irregularities from
differential ice core melting (Figure A.5). This suggests that thaw-related processes likely
contribute minimally to geologic scatter in the Itsi Range.

38

Subglacially entrained debris is of equal or greater importance than supraglacially
entrained debris in low-relief glaciated valleys, significant because a fundamental assumption of
TCND is that the rocks are subglacially, not supraglacially, i.e., rockfall, derived (Barr and
Lovell, 2014). Rockfall derived supraglacially transported boulders bearing inherited nuclides
and deposited as a moraine would represent maximum ages that also do not reflect the true
depositional age (Hawkins, 2023). All lines of evidence considered, recent literature suggests
post-depositional modification plays a less significant role for cirque glaciers than boulders
derived from plucking of pre-exposed boulders from a glacier valley and boulders reincorporated
from adjacent rockfall (Goehring et al., 2022; Tomkins et al., 2021).
2.6.2 Inner moraine
Boulders sampled from the inner moraine yield a median exposure age of 0.20 ka ± 0.05
(Table 2.7). The dating uncertainty of glacier advances generally increases with age, so recent
advances (e.g., LIA) are more confidently dated than older (e.g., LGM) ones (Kirkbride and
Winkler, 2012). For the inner moraine, I assert that the large X2 value of 12.73 does not indicate
the presence of significant geological scatter, rather, it reflects the small sample size and narrow
range of uncertainty associated with the extremely young ages. Exceptionally young surface
exposure ages have narrow uncertainty windows. This small uncertainty, coupled with the
heightened landform instability of young moraines (Tomkins et al., 2021), makes the ages less
likely to overlap. The MAD value means that, on average, the exposure ages for the inner
moraine deviate from the median by 87 years. This represents a close clustering of ages around
the median; since the inner moraine IQR is also low, the MAD reflects a narrow time range of
glacial activity.

39

Conversely, strict interpretation of the reduced chi-squared statistic of 12.73 implies that
geomorphological factors beyond random measurement errors affect the surface exposure ages.
For young moraine landforms (<1 ka), steep ice-proximal slopes reflect instability of the
inchoate landform, which makes sampling moraine crests a poor strategy if making efforts to
minimize boulder instability (Tomkins et al., 2021). An improved strategy of future sampling
recommendations for presumably young, latest Holocene landforms should avoid solely
sampling from moraine crests and instead sample from alternative intra-landform locations of
ice-proximal and ice-distal slopes.
2.6.3 Outer moraine
Boulders sampled from the outer moraine yield a median exposure age of 0.96 ka ± 0.36
(Table 2.7). This moraine displays obvious signs of geological scatter, returning a X2 of 133.44,
further supported by a MAD of 362 years. The anomalously old ages probably reflect
accumulations of inherited nuclides. Melting permafrost in ice-cored moraines likely did not
cause thaw slumps or contribute to the observed scatter, as boulder ages skew older than
expected for Holocene maximum extents. Moraine boulders derived from plucking of preexposed boulders from a glacier valley and boulders reincorporated from adjacent rockfall may
lead to inheritance in a cosmogenic nuclide dated sample. Regional observations of seismically
induced large-scale rockslides further support mechanisms of inheritance rather than thaw slump
induced post-depositional disturbance (Smith and Bennett, 2004).
Observations in the field (Figure 2.8) and analysis of the 2004 DEM (Figure 2.4) show an
elongated ridge above the inner moraine, proximal to the ice-margin, signifying exposed rock
repeatedly abraded during periods of less extensive glacier limits. This rock could have been
subsequently overridden and redeposited at the existing terminal moraine (Balco, 2020a). If we
40

assume that the exposure ages are accurate which are a consequence of the actual climatic
processes that led to deposition of the moraine, it must be considered that the moraine consists of
an amalgamation of prior overridden advances; a representation of multiple neoglacial advances
and retreats over the past two millennia.

Figure 2.8: An oblique angle of the Itsi Glacier forefield highlighting the inner and outer
moraines as well as an elongated ridge of glacially abraded rock proximal to the ice margin and
parallel to the moraine axes.
2.6.4 Erratics
Sampled erratic boulders yield a median exposure age of 11.61 ka ± 0.54 (Table 2.7)
denoting precise accordance of CIS melt with the end of the YD stadial recorded globally
(Broecker et al., 2010; Menounos et al., 2017). This suggests that as the CIS deteriorated, it
transitioned into a large valley glacier evidenced by the deposition of a large moraine far
downvalley of the Itsi Glacier outer moraine. The valley glacier subsequently thinned and melted
away, depositing the dated erratics immediately beyond the Itsi Glacier outer moraine. This
indicates that the CIS was no more extensive than the erratics’ position at this time. Both erratics

41

fail the reduced chi-squared test albeit with a low chi-squared value of 3.68. The very high MAD
value of 1073.5 years highlights extreme deviation from the median, but the small sample size (n
= 2) limits ability for strict interpretation and indicates poor statistical robustness. Considering
this, the boulders returned ages consistent with the regionally interpreted timeline of CIS decay,
implying sample 228B10 was negligibly affected by the steep extraction angle (Menounos et al.,
2017). This suggests that, at least at this site, the extraction angle may be less significant in
influencing exposure ages than previously thought.
2.6.5 Geodetic mass balance and termini retreat
The large standard deviations associated with glacier terminus mean retreat rates reflect
both the morphological diversity among glacier termini and the method used for their calculation
(Urbanski, 2018). ERA5-Land data for the Itsi Range from 1950 to 2023 helps contextualize the
significant increase in retreat rate and mass loss since deposition of the Itsi Glacier outer moraine
(Tables 2.5 and 2.6). Glaciers in interior continental regions at high altitudes, such as the Itsi
Glacier, primarily respond to fluctuations in summer temperatures during the melt season (Roe
and O’Neal, 2009). Boreal summer and winter trends reveal a steady increase in mean
temperature and precipitation along with a decrease in summer snow cover since 1950 (Figure
2.9). Although an increase in mean temperature is the most direct contributor to termini retreat
and loss of glaciated area among Itsi Range glaciers, increasing summer precipitation on
diminishing snow cover can further enhance glacier surface melting. MAT anomalies for 2023
provide additional evidence of this accelerated warming trend (Figure 2.10).

42

Figure 2.9: ERA5-Land values for the Itsi Range (62.9°, -130.2°) plotted to represent boreal
summer melt (June, July, August (JJA)) and winter accumulation (December, January, February
43

(DJF)) seasonal means. Values are plotted as five-year moving averages with associated
trendlines.

Figure 2.10: ERA5-Land 2023 MAT anomalies for the YT using reference period 1950-1980.
Itsi Range is highlighted in green.
The Itsi Range straddles the Upper Pelly and Macmillan drainage basins within the
greater Pelly catchment. Mean monthly discharge observations from the Ross River station,
situated downstream of the Itsi Range in the Upper Pelly basin, corroborates with modelled
ERA5 temperature and snow cover trends. Discharge trendlines illustrate a shift to an earlier
melt season peaking in May, a departure from historically consistent observations of peak
discharge values in June (Figure 2.11) (Smith and Bennett, 2004). This shift results in an earlier
disappearance of snow cover on Itsi Range glaciers and consequently earlier exposure of bare
ice. Exposed glacier ice accumulates light-absorbing impurities (e.g., algae, dust, supraglacial

44

meltwater, etc.) which reduces overall albedo, and, in combination with a prolongation of melt
season as indicated by discharge trends, amplifies mass loss (Naegeli et al., 2019).

Figure 2.11: Observed mean discharge values at Ross River (61.99, -132.41; Station ID:
09BA001) within the Upper Pelly sub-sub-drainage area for May and June. Extracted from the
Environment and Climate Change Canada Historical Hydrometric Data web site
(https://wateroffice.ec.gc.ca/mainmenu/historical_data_index_e.html).
I compared the geodetic mass balance records obtained for the Itsi Range with those
published by Hugonnet et al. (2021), utilizing the same glaciers to assess the consistency of
volumetric and mass change values (Table 2.8). The results from this study show a volumetric
change of -0.091 km3 ± 0.034 over 18 years, with a corresponding mass change rate of -0.382 m
w.e. yr⁻1 ± 0.136 across an analysis area of 11.180 km2. In comparison, Hugonnet et al. (2021)
reported a volumetric change of -0.180 km3 ± 0.331 over 15 years, with a mass change rate of 0.692 m w.e. yr⁻1 ± 0.649 for a larger analysis area of 15.022 km2.
Table 2.8: Mass budget of Itsi Range glaciers compared to Hugonnet et al. (2021).
Source
This study
Hugonnet et
al. (2021)
a

years
19

Analysis
ΔV
area
(km2)
(km3)a
11.180 -0.091 ± 0.034

(m w.e. yr-1)a
-0.382 ± 0.136

16

15.022

-0.649 ± 0.609

Start

End

Elapsed

year
2004

year
2022

2004

2019

-0.180 ± 0.331

Mass change

Uncertainty estimations for 𝐸Δ𝑉 and σΔ𝑀,𝑟𝑑𝑛 match those described in section 2.4.2.

45

A portion of the differences in the volumetric and mass change values may be attributed
to the sources used for analysis. Hugonnet et al. (2021) use the Randolph Glacier Inventory (RGI
6.0) in tandem with Advanced Spaceborne Thermal Emission and Reflection Radiometer
(ASTER) DEMs. In contrast, I created dynamic glacier outlines for each year and processed
higher resolution (8.11 m) DEMs sourced from aerial photographs and ArcticDEM. The RGI
delineations are inherently less accurate due to the automated nature of its creation and static
area over the assessed time period. Additionally, the 30 m resolution of ASTER DEMs heightens
potential for error. Despite these methodological differences, both studies corroborate the overall
trend of accelerated glacier mass loss in the Itsi Range and suggest that this study’s geodetic
mass balance estimations are conservative. These trends in mass loss suggest that, due to the
lagged glacier response to climate fluctuations, retreat will continue even if the increase of
temperature anomalies ceased, as glaciers in the Itsi Range remain in disequilibrium with the
warming climate (Marshall et al., 2011).
2.7 Conclusion
I used TCND techniques on a pair of nested moraines to develop a latest Holocene
chronology of the Itsi Glacier among the Itsi Range in the YT. Erratics from this site-specific
chronology suggests that the decaying CIS exposed Itsi Range cirque glaciers correlative to the
end of the YD, meaning that cirque glacier advances at this time were less extensive than their
LIA extents. The most extensive cirque glacier advances occurred during the latest Holocene,
recorded as two distinct pulses and further evidenced by the existence of several nested moraine
pairs throughout the massif. Inner moraine exposure age of 0.20 ka ± 0.05 is minimally affected
by geologic scatter with uncertainty a reflection of moraine crest instability due to the young age
of the landform, limited sample size, and narrow ranges of age uncertainty among individual
46

samples. Conversely, the outer moraine age of 0.96 ka ± 0.36 exemplifies the issues of geologic
scatter commonly reported in the northeastern Canadian Cordillera. The lack of further evidence
allows for the following conclusions to be drawn about the outer moraine: 1) the oldest boulder
age of 2.09 ka ± 0.15 is correct and corresponds to the time of moraine construction, with
younger boulder ages reflecting subsequent moraine stabilization, 2) the moraine reflects the
thaw slump of an ice-cored moraine causing erroneously young ages, 3) the moraine reflects
processes of inheritance due to plucking of pre-exposed boulders from the glacier forefield
and/or from adjacent rockfall causing erroneously old ages, 4) all boulder ages are true and
accurately reflect regional climatic processes, indicating an accretion of overridden glacier
advances.
GTT derived rates of glacier retreat and geodetic mass balance results for Itsi Range
glaciers point to an escalating trend of retreat in recent decades attributed to a consistent increase
in MAT and decrease in summer snow cover since 1950. Discharge observations downstream
from the Itsi Range since 1960 capture evidence of a shift to an earlier melt season which will
propagate trends of glacier mass loss in this region. ERA5-Land temperature and snow cover
trends further substantiate the effects of rapid warming that began in the 20th century.

3. Comparison with western North America moraine records
3.1 Introduction
The evidence of Holocene glacier fluctuations within the northeastern Canadian
Cordillera remains sparse. Chapter 3 seeks to compare the Itsi Glacier 10Be chronology with
available regional records of moraine and erratic ages. In addition, I provide complementary preexisting paleoclimate chronometer and proxy records from radiocarbon, dendrochronology, and
47

lake sediment cores to further improve understanding of how these glacial landforms fluctuated
throughout the latest Holocene.
Both the frequency and magnitude of climate drivers and the rate and duration of glacier
advance contribute to the moraine record glaciers emplace over time (Kirkbride and Winkler,
2012). A theoretical experiment conducted by Kirkbride and Winkler (2012) found that glaciers
with short response times to fluctuations of climate are more sensitive and record more advances
in the moraine record than glaciers with long response times, wherein the process of obliterative
overlap begins to censor the moraine record. If multiple glaciers have a similar age for a given
moraine, it signifies an advance likely occurred followed by a long period of retreat, decreasing
the likelihood of censoring by obliterative overlap (Kirkbride and Winkler, 2012). This reasoning
provides the paradigm for comparing regional glacial records to draw paleoclimatic conclusions,
as a lower likelihood of erosion censoring reduces opportunities for false correlations among
geographically discrete glacier advances. Obliterative overlap, in the form of late Holocene
glacier expansion, is the most likely reason for the lack of preservation of pre-Little Ice Age
(LIA) moraines at the Itsi Glacier site.
3.2 Methodology
I primarily sourced terrestrial in situ cosmogenic nuclide (TCN) dated moraine boulder
ages from the Informal Cosmogenic-Nuclide Exposure-age Database (ICE-D): Alpine project
(Balco, 2020b). As few TCND studies have been carried out in the northeastern Canadian
Cordillera, this chapter incorporates all available late Holocene 10Be dates from moraines in
Alaska (AK), British Columbia (BC), Yukon Territory (YT), and Northwest Territories (NWT). I
first recalculated 10Be samples to modern standards matching those described in chapter 2 to
ensure consistency across regional comparisons.
48

I assessed all relevant 10Be dated moraine positions in Western North America (WNA)
and assigned them one of four positions as defined in Table 3.1. Rather than relying solely on
10

Be-dated moraine ages, morphostratigraphic position was the primary criterion for assigning

end groups, with moraine age the secondary criterion. This approach minimizes the impact of
geologic scatter on moraine age estimates and allows for the identification of potential patterns in
scatter associated with morphostratigraphic position. I normalized reported ages to ka BP and
present them as median ± interquartile range (IQR). To assess the presence of geologic scatter, I
first applied the KS test to determine whether a parametric or non-parametric statistical approach
was more appropriate. For parametric groups I report the reduced chi-squared statistic, whereas
the MAD is provided for non-parametric groups.
Table 3.1: Description for each morphostratigraphic position used for WNA comparisons.

a

Position

Agea

Sample locationsb

End
Group 1

LIA (~0.90 – 0.06 ka)

AK, BC, NWT,
YT

End
Group 2

LIA (~0.90 – 0.06 ka)

AK, NWT, YT

End
Group 3

Pre-LIA (~4.40 – 1.20 ka)

AK, BC

Erratics

CIS Decay (~12.7 ka)

BC, NWT, YT

Description
Singular innermost recessional
or lateral moraine proximal to
ice-margin.
Any recessional or lateral
moraines distal to End Group 1.
Any recessional or lateral
moraines distal to End Group 1
and/or 2.
Erratic boulders.

End Group 1 – 3 ages estimated using WNA Neoglacial advances discussed in chapter 3.3.1.
AK = Alaska, BC = British Columbia, NWT = Northwest Territories, and YT = Yukon Territory.

b

In addition to 10Be surface exposure ages, I grouped and assessed proxies derived by
methods of radiocarbon dating, dendrochronology, and glaciolacustrine sedimentation records.
These proxies, synthesized from numerous publications, correspond to intervals of glacial
advance reported in the literature (Barclay et al., 2009; Menounos et al., 2009; Mood and Smith,
2015; Tomkins et al., 2008).

49

Radiocarbon dating techniques provide ages for when a carbon-based organism has died,
mainly as subfossil wood (Braumann et al., 2020; Mood and Smith, 2015). Radiocarbon dating’s
relevance in reconstructing glacier chronologies relies on the assumption that advancing glaciers
may kill trees growing within glacial valleys, preserving them as glaciogenic deposits. During
subsequent glacier retreat, overridden wood becomes exposed. Once exposed, the sampled wood
can provide a radiocarbon-derived kill date assigned to a glacier advance (Ryder and Thomson,
1986).
Dendrochronology, utilizes the growth of tree-rings, useful for resolving high temporal
resolution glacier histories during the LIA (Jacoby and Cook, 1981; Osborn et al., 2001;
Tomkins et al., 2008). The dendrochronology dating method works by analyzing the patterns of
tree ring widths and ring densities, which vary annually due to a variety of environmental
conditions. Tree rings record an aggregation of multiple environmental changes; therefore,
depending on the application, e.g., correlating with summer temperature records, the desired
signal must first be extracted from the noise of irrelevant contributing factors (Cook and
Kairiukstis, 1990).
Glacially distal lacustrine sediments present a high resolution proxy record that records
influxes of sediment that reflect regional environmental changes helpful in reconstructing
regional paleoclimate and related glacial histories (Leonard, 1997; Nesje et al., 2001; Tomkins et
al., 2008). Grain size variations in glaciolacustrine sediment records reflect the activity of wetbased, temperate glaciers (Nesje et al., 2001). Organic content is estimated using percent weight
loss-on-ignition (LOI), which measures sample weight before and after high-temperature
exposure. Increased erosion and input of clastic sediments, indicating heightened glacial activity,
can be inferred from increases in magnetic susceptibility. Together, LOI and magnetic
50

susceptibility data enable continuous intercalation of Holocene glacier fluctuations (Maurer et
al., 2012; Nesje et al., 2001). Although radiocarbon dating and dendrochronology primarily
provide limiting ages, when combined with glaciolacustrine sedimentation records and/or 10Be
dating, they enhance the precision of glacier fluctuation chronologies by better constraining the
time series of glacier advances (Hawkins et al., 2021; Osborn et al., 2007).
For climatic context, I review the timing of alpine glacier advances among WNA from
the early Neoglacial to the present, drawn from these published sources. I compiled these
advances to represent minimum and maximum ages of advance (e.g., for a bracketed age of 1.9 –
1.4 ka, 1.9 ka would be the maximum-limiting age). Periods of glacial advance over the past
millennia are presented visually using Python (v. 3.11.7).
3.3 Discussion of moraine records
3.3.1 Review of western North America Neoglacial advances
The Northern Hemisphere (NH) experienced generally cold, glacially extensive period
following the Holocene Thermal Maximum (HTM); however, there still exists a significant
degree of intrinsic climate variability (Clague et al., 2009; Renssen et al., 2009). Radiocarbon
ages record two distinct early-middle Neoglacial advances between 4.40 – 3.97 ka and 3.70 –
3.40 ka in the western Canadian Cordillera (Menounos et al., 2009; Mood and Smith, 2015).
Generally, advances were slightly less extensive during the middle-to-late Neoglacial than during
the LIA. Western Canadian Cordillera and Alaskan (AK) glaciers record three distinct advances
during this period evidenced by dendrochronology and radiocarbon kill dates at 3.54 – 2.50 ka,
2.20 – 2.00 ka, and 1.71 – 1.20 ka. Literature frequently refer to the advance between 3.54 – 2.50
ka as the contemporaneous Tiedemann-Peyto Advance of the BC Coast and Rocky Mountains

51

(Luckman et al., 1993; Menounos et al., 2009; Ryder and Thomson, 1986). The TiedemannPeyto Advance is complex, more accurately characterized by multiple periods of glacier advance
and retreat (Menounos et al., 2009). In AK, Barclay et al. (2009) reported radiocarbon evidence
for an advance at 2.20 – 2.00 ka. Dendrochronology and radiocarbon evidence exists for a third
advance interval of similar complexity to the Tiedemann-Peyto Advance, often referred to as the
First Millennium Advance, throughout the western Canadian Cordillera and AK occurring during
the interval 1.71 – 1.20 ka (Barclay et al., 2009; Menounos et al., 2009; Mood and Smith, 2015;
Reyes et al., 2006). St. Elias Mountain valley glaciers advanced from 3.30 – 2.40 ka, coeval with
the Tiedemann-Peyto Advance; this contrasts with younger ages reported within this region
during the First Millennium Advance from 1.25 – 1.05 ka (Menounos et al., 2009). The large
interval of time ascribed to the First Millennium Advance suggests complexity of glacial
behavior akin to the Tiedemann-Peyto Advance (Menounos et al., 2009).
Neoglacial glacier advances occurred globally during the LIA, but there is great spatial
variability in the timing of these advances at the centennial scale (Kaplan et al., 2016; Schaefer et
al., 2009; Solomina et al., 2015). Many Canadian Cordillera glaciers attained their Holocene
glacial limits during the mid-late 19th century (Menounos et al., 2013a, 2009; Solomina et al.,
2015). Dendrochronology, radiocarbon, and glaciolacustrine evidence identifies three phases of
advance throughout WNA, all generally referred to as the LIA, from 0.80 – 0.40 ka, 0.50 – 0.24
ka, and 0.40 – 0.06 ka (Barclay et al., 2009; Hawkins et al., 2021; Luckman, 2000; Menounos et
al., 2009; Mood and Smith, 2015; Solomina et al., 2015; Tomkins et al., 2008). The compilation
of major WNA advances are presented in Table A.1

52

3.3.2 Radiocarbon, dendrochronology, and glaciolacustrine records
When a glacier enters equilibrium with the predominating climate, emplacement of
terminal and lateral moraines occurs (Solomina et al., 2015). This means that TCND derived
moraine ages do not directly indicate when the glacier started advancing; rather, it signifies the
culmination of advance, i.e., onset of retreat. Itsi Glacier moraine ages correlate variably with
glacial advances across WNA (Figure 3.1).

Figure 3.1: Regional compilation of glacial advance intervals from 4.0 ka to present based on
regional reviews using evidence from radiocarbon dating, dendrochronology, and
glaciolacustrine sedimentation records (Barclay et al., 2009; Hawkins et al., 2021; Menounos et
al., 2009; Mood and Smith, 2015; Tomkins et al., 2008). Grey box widths are set by duration of
advance interval with ages calibrated to ka BP 1950. Dashed lines indicate approximate age
intervals when discrete ages were not provided. Orange and blue vertical bars are Itsi Glacier
median moraine ages, signifying moraine stabilization and glacier retreat; bar widths are set by
the IQR of 10Be dated moraine ages. Pink vertical bars denote ages of individual Itsi Glacier
outer moraine boulders which comprise the median outer moraine age; bar widths are set by the
propagated error of each boulder.
The outer moraine dates to approximately 250 years after the First Millennium Advance,
1.71 – 1.20 ka (Figure 3.1). The large IQR and scattered exposure ages suggest multiple
possibilities for moraine genesis as discussed in chapter 2. One explanation is that this moraine
includes boulders from smaller moraines that were destroyed when the glacier attained its
53

maximum Holocene extent. Extending the glacial advance record and plotting individual outer
moraine boulder ages reveals the complexity of late Holocene advances in this region. The visual
correlation in Figure 3.1 establishes a plausible relationship between glacial advances and outer
moraine boulders. Notably, the onset of three advances precedes outer moraine boulder ages
when accounting for error. In Alaska and western Canada, the record of middle-late Neoglacial
regional advances at 2.20 – 1.80 ka and 1.71 – 1.20 ka align well with the outer moraine boulder
ages of 2.09 ka ± 0.15 and 1.16 ka ± 0.08, respectively (Barclay et al., 2009; Hawkins et al.,
2021; Menounos et al., 2009; Mood and Smith, 2015). Additionally, the latest Neoglacial
regional advance from 0.90 – 0.40 ka coincides with the outer moraine boulders dated to 0.43 ka
± 0.04 and 0.76 ka ± 0.06 (Barclay et al., 2009; Hawkins et al., 2021; Menounos et al., 2009;
Mood and Smith, 2015; Solomina et al., 2015). Glacial lake sediments assessed by Tomkins et
al. (2008) indicate an advance among the southern Selwyn Mountains that ended at 0.50 ka, but
the record does not extend prior to 0.56 ka, leaving timing for the onset of this advance
uncertain. The advance likely correlates with other regional records beginning around 0.90 –
0.77 ka (Barclay et al., 2009; Hawkins et al., 2021; Mood and Smith, 2015).
The Holocene maximum extent of the Itsi Glacier occurred at 0.96 ka ± 0.36. In contrast,
Tomkins et al. (2008) cite Holocene maxima in the southern Selwyn Mountains at 0.17 – 0.06
ka. Menounos et al. (2009) and Mood and Smith (2015) similarly interpret later maximum
Holocene extents during the 19th century. Hawkins et al. (2023) identify an earlier Holocene
maximum at 0.39 ka. Both the Itsi Glacier and the sites sampled by Hawkins et al. (2023) in the
Selwyn and Mackenzie Mountains display nested moraine formations beyond alpine cirques.
Sample 228B05 from the Itsi Glacier end group 2 moraine provides an age of 0.43 ka ± 0.04,
plausibly making this the true minimum-limiting age for LIA maximum extents if accepting the
54

interpretation that the moraine represents an amalgamation of multiple glacial advances; or, if
interpreting the other end group 2 samples as containing inherited nuclides. Evidence for a circa
0.40 ka Holocene maximum is further supported by Barclay et al. (2009), who report southern
coastal AK glaciers attained Holocene maximum extents during both middle and late LIA
advance phases circa 0.41 – 0.24 ka and 0.14 – 0.07 ka. Given the limited glacier chronology
records, this implies that in the northeastern Canadian Cordillera, LIA maximum extents were
reached during the 0.90 – 0.40 ka advance, centuries earlier than in the western and southwestern
Canadian Cordillera.
The Itsi Glacier inner moraine was emplaced at 0.20 ka ± 0.05 and is interpreted to
represent the culmination of a singular regional advance from 0.40 – 0.06 ka. Specifically, Itsi
Glacier inner moraine ages correspond closely to an advance in the southern Selwyn Mountains
prior to onset of glacial recession from 0.22 – 0.17 ka (Tomkins et al., 2008). Coeval intervals of
advance are observed in the Coast Mountains of BC at 0.24 – 0.07 ka and the Sierra Nevada
Mountains of the American Cordillera at 0.25 – 0.17 ka (Mood and Smith, 2015; Solomina et al.,
2015). Varved lake sediments in the southern Selwyn Mountains record a later glacial recession
starting at 0.06 ka, suggesting inner moraine surface exposure ages represent maximum limiting
ages for the culmination of the glacial advance (Tomkins et al., 2008).
3.3.3

10Be moraine ages

Hawkins et al. (2023) reported consistent LIA exposure ages circa 0.30 and 0.10 ka from
outer and inner moraine pairs in the northeastern Canadian Cordillera, along with outliers
ranging from 0.99 ka to 4.31 ka. Similar observations of scatter have been reported for late
Pleistocene moraines along the northwestern margin of the Cordilleran Ice Sheet (CIS) and
among mountain glaciers of the Grey Hunter massif in the YT, making clear interpretations of
55

these ages difficult (Goehring et al., 2022; Stroeven et al., 2014). A critical question in the field
of TCND is whether these outlier ages represent true landform ages or geological scatter. The
MAD is used for discussion of geologic scatter as all end groups reject the null hypothesis of
normality as posed by the KS test (Table 3.2).

Figure 3.2: Regional compilation of WNA late Holocene 10Be boulder ages, including those
collected for this thesis, sampled from moraines and grouped by morphostratigraphic position
(Dortch et al., 2010; Hawkins et al., 2023, 2021; Howley, 2008; Menounos et al., 2013a;
Pendleton et al., 2017; Young et al., 2009). Thin grey curves show ages of individual boulders
rendered as Gaussian distributions. Thick curves represent the sum of each Gaussian distribution
for the entire population of 10Be moraine group ages, organised by colour. Itsi Glacier moraine
groups are inclusive in the regional curves and also highlighted separately within the summed
Gaussian distributions. Corresponding box-and-whisker plots denote the median ± IQR,
minimum, and maximum surface exposure ages of each moraine group.
Itsi Glacier end group 1 ages accord well with the regional 10Be dated moraine record
(Figure 3.2 and A.1). Although we did not sample lateral moraines for this study, I included them
in end group 1 to provide better context for the complex timing of glacial standstills/retreats in
WNA during the last millennium. 10Be moraine ages throughout the Canadian Cordillera (n =
21) and AK (n = 4) provide median end group 1 ages of 0.13 ka ± 0.1 (n = 25). The KS test

56

yielded a p-value of less than 0.05, indicating End Group 1 is a non-parametric distribution of
ages (Table 3.2). The MAD of 67 years reflects very low variability among the exposure ages in
this group, suggesting that the ages in End Group 1 likely correspond to a specific climatic event.
End group 1 moraines sampled at high latitudes (62.1° – 63.6°) across AK, NWT, and
YT date to 0.30 ka ± 0.10 (n = 11) (Dortch et al., 2010; Hawkins et al., 2023). At southern
latitudes (50.9°), lateral moraines in the BC Coast Mountains record end group 1 emplacement at
0.09 ka ± 0.03 (n = 14) (Hawkins et al., 2021). End group 1 moraines show less geological
scatter than end group 2, the latter returning an MAD of 201 years (Table 3.2).
Table 3.2: 10Be dated WNA moraines grouped by morphostratigraphic position with associated
statistics (Dortch et al., 2010; Hawkins et al., 2023, 2021; Howley, 2008; Menounos et al.,
2013a; Pendleton et al., 2017; Young et al., 2009).
Position
End group 1
End group 2
End group 3

n
25
36
42

Median (ka) IQR (ka)
0.13
0.10
0.46
0.35
2.05
0.69

X2
10.82
129.18
321.50

MAD
67
201
787

D-stata
0.283
0.303
0.436

p-valueb
0.029
0.002
0.000

a

KS test D-statistic: the maximum difference between the empirical CDF of the sample and the CDF of the normal
distribution (Massey Jr, 1951).
b
KS test p-value: the probability that the sample comes from a normal distribution. A value less than 0.05 indicates
a rejection of the null hypothesis, indicating data are not normally distributed (Massey Jr, 1951).

Itsi Glacier end group 2 ages fall within the 2σ uncertainty represented by the Gaussian
distributions, skewing older than the regional record (Figure 3.2). Hawkins et al. (2023) dated
eight outer moraines from paired LIA moraines among the Mackenzie and Selwyn mountain
ranges. Out of 23 moraine group 2 samples, eight returned anomalously old 10Be ages,
purportedly due to inheritance, ranging from 0.99 ka ± 0.10 to 4.31 ka ± 0.29. This pattern is
consistent among end group 2 moraines in the Mackenzie and Selwyn Mountains, with all
anomalously old ages in the region possessing a combined median of 2.29 ka ± 0.83 (n = 10).
These ages closely resemble AK end group 3 10Be ages at northern latitudes (63° – 68°), which
have a combined median of 2.33 ka ± 0.63 (n = 20) (Pendleton et al., 2017; Young et al., 2009),
57

and an end group 3 moraine at the Tiedemann Glacier in BC dated to 2.60 ± 0.28 (n = 6)
(Menounos et al., 2013a).
End group 3 gaussian distributions illustrate significant overlap with end group 2, largely
due to instances of geologic scatter among end group 2 moraines. The MAD values of 201 and
787 years for end groups 2 and 3 respectively. The exposure ages in these groups do not
represent a coherent climatic event and may reflect either a broader range of geological
processes, such as inheritance or post-deposition disturbance, or multiple climatic events
affecting the moraine formation. Alpine glaciers in northern and interior AK, as well as glaciers
among the Wrangell and Saint Elias mountains, contain pre-LIA moraines rarely preserved
within WNA due to obliterative overlap achieved during LIA maximum glacial extents (Barclay
et al., 2009; Denton and Karlén, 1977; Pendleton et al., 2017). In the Brooks Range of AK,
moraines overridden by pre-LIA advances record maximum Holocene extents between 3.5 – 2.6
ka (Pendleton et al., 2017). Subsequent advances through to the LIA are evident, but cirque
glaciers in the Brooks Range deposit topographically high moraines, making it difficult for the
short, debris-rich glaciers to override them (Pendleton et al., 2017). The lack of LIA moraines
may indicate morphological constraints inherent to the glaciers themselves rather than climatic
trends, i.e., the lack of LIA maximum extents does not reflect a warming climate (Denton and
Karlén, 1977; Pendleton et al., 2017). These pre-LIA Brooks Range moraines are coeval to the
complex Tiedemann-Peyto Advance, suggesting a relation to multiple minor pulses of advance
between 3.54 – 2.50 ka (Menounos et al., 2013a).
End group 2 moraines in the northeastern Canadian Cordillera represent the most
extensive glacial limits throughout the Holocene, emplaced in response to the 0.90 – 0.40 ka
regional advance. These moraines likely contain remnants of early Neoglacial advances reflected
58

by outlier 10Be ages that have been merged. Specifically, the Itsi Glacier and other end group 2
moraines potentially created a complex, amalgamated landform due to the inability of small
alpine glaciers to override topographically high barriers formed by early Neoglacial terminal
moraines. This amalgamation mechanism suggests that these glaciers, when advancing during
the Holocene maximum extent circa 0.40 ka, encountered pre-existing moraines from earlier
neoglacial phases of similar magnitude (Figure 3.4). Instead of completely overriding and
resetting the 10Be boulder inventory, these glaciers added new material to the existing moraine.
This process preserved older boulders with inherited 10Be ages, resulting in the consistent reports
of geological scatter.

Figure 3.3: Regional compilation of WNA late Holocene 10Be moraine ages, including those
collected for this thesis, grouped by morphostratigraphic position (Dortch et al., 2010; Hawkins
et al., 2023, 2021; Howley, 2008; Menounos et al., 2013a; Pendleton et al., 2017; Young et al.,

59

2009). Left: map illustrating location of sampled moraines where size of marker is proportional
to IQR. Right: scatter plot showing median moraine ages in relation to IQR.
In contrast, end group 1 moraines, culminating during the 0.40 – 0.06 ka regional
advance, did so on relatively “clean” surfaces devoid of significant prior moraine deposits.
Although the sample size of 10Be dated moraines in this region is small (n = 26), a discernible
pattern emerges (Figure 3.3 and A.2). End group 2 moraines generally exhibit a consistently
higher IQR compared to both end group 1 and, in most cases, end group 3. Moreover, the highest
IQR values are centered along the YT/NWT border. The lack of inherited boulders from earlier
advances resulted in more consistent 10Be ages and minimal scatter. The morphological
constraints and subsequent amalgamation of moraines in the northeastern Canadian Cordillera
provide a plausible explanation for the significant age scatter observed in end group 2 moraines.

60

Figure 3.4: Conceptual model of moraine amalgamation adapted from Bierman et al. (2021). An
amalgamated moraine is an unstratified mixture of clasts, some of which have experienced
periods of glacier transport and deposition one or several times before final transport and
deposition on the contemporary moraine. At time step 𝑡1 , (i) buried clasts are shielded from
cosmic rays, (ii) glacially transported, and (iii) emplaced as a terminal moraine following
deglaciation. At 𝑡2 , (iv) buried boulders are, (v) glacially transported, and (vi) emplaced as a
recessional moraine proximal to the ice margin. At 𝑡3 , the glacier (vii) overrides the prior
moraine deposit from the second glaciation and (viii) emplaces an unstratified mixture of
previously deposited clasts as an amalgamated terminal moraine. The resulting exposure age
distribution (ix) reflects multiple records of deglaciation.
3.3.4

10Be erratic ages

I include additional ages from erratics and bedrock boulders to complement cirque and
valley moraine ages for constraining the timeline of CIS decay (Menounos et al., 2017). Itsi
Glacier erratics indicate CIS melt exposed alpine glaciers circa 11.61 ka ± 0.54 (n = 2) among
the Itsi Range. This accords with the larger dataset of Canadian Cordillera erratic samples 12.13
ka ± 1.69 (n = 26) and resembles observations by Hawkins et al. (2023) in the northeastern
Canadian Cordillera who recorded erratic deposits ranging from 10.94 ka ± 0.31 to 12.27 ka ±
0.22 (n = 4) (Figure 3.5). Together, these erratics date to 11.56 ka ± 1.55 (Table 3.2). Alpine
valley and cirque glaciers advanced coeval with and preceding the Younger Dryas (YD) interval
throughout the Canadian Rockies (Hawkins et al., 2023; Luckman and Osborn, 1979; Menounos
et al., 2009) which suggests cirque glacier advances correlative to the YD among the Itsi Range
were less extensive than their LIA extents. Given that the suite of WNA erratic samples reflect
the high spatiotemporal variability of CIS decay over millennia, it does not correspond to a
single climatic event; thus, no attempt is made to address scatter in this context.

61

Figure 3.5: Regional compilation of WNA 10Be erratic boulder ages, including those collected
for this thesis, from BC and YT (Hawkins et al., 2023; Menounos et al., 2017; Stroeven et al.,
2010). Thin grey curves show ages of individual boulders rendered as Gaussian distributions.
The thick pink curves represent the sum of each Gaussian distribution for the entire population of
10
Be erratic ages. Itsi Glacier erratics are highlighted separately within the summed Gaussian
distributions. A corresponding box-and-whisker plot denotes the median ± IQR, minimum, and
maximum surface exposure ages for all erratics.
3.4 Limitations
Different geochronological methods have their own formation processes and variable
sensitivity to their environment (Tomkins et al., 2008). I compiled intervals of advance derived
from radiocarbon, dendrochronology, and glaciolacustrine sedimentation methods (Barclay et
al., 2009; Hawkins et al., 2021; Menounos et al., 2009; Mood and Smith, 2015; Tomkins et al.,
2008). As they do not record identical climate signals, evaluating these different methods for
timing of glacier advance and moraine emplacement requires consideration of their limitations.
Glaciolacustrine sedimentation records are a proxy for glacial meltwater production,
which is dictated by regional climate and glacier dynamics (Tomkins et al., 2008). Glaciers have
a lagged response to changes in climate. The magnitude of this lag can vary by tens to hundreds
of years between mountain glaciers. Non-climatic factors affecting proglacial sedimentary record
62

formation, like sediment supply, can obscure the glaciogenic clastic signal and consequent
interpretation (Menounos et al., 2006; Tomkins et al., 2008). Although sedimentation rates
generally increase with glacial ice extent, maximum ice extents do not always align with the
highest rates of sedimentation (Leonard, 1997). The relationship between clastic sedimentation
rates and glacier extent is more consistent on centennial-to-millennial timescales (Leonard,
1997). On shorter, decadal timescales, the highest sedimentation rates often occur during
transitional periods, such as during glacier advances before a maximum or recessions afterward,
possibly obscuring the signal of glacier advance (Leonard, 1997).
Dendrochronology records reflect immediate changes in climate. Tree-ring patterns can
also record non-climatic factors like fire, disease, and tree physiology, making it difficult to
disaggregate the desired signal (Cook and Kairiukstis, 1990; Tomkins et al., 2008).
Radiocarbon dating techniques have two main sources of uncertainty: 1) high uncertainty
of radiocarbon ages relative to absolute age within the last millennium and 2) temporal and
spatial uncertainties. Temporally, radiocarbon dates can only bracket a given glacier advance,
and spatially, they rarely provide precise glacier positions unless dated materials are found in situ
and still in growth position (Menounos et al., 2009; Mood and Smith, 2015). Climatic signals can
still be ascertained among a suite of different chronometer and proxy measurements if they are
all similarly influenced by the same climatic factors of sufficient magnitude (Braumann et al.,
2020; Tomkins et al., 2008).
I synthesized morphologically relevant 10Be moraine ages because Balco (2020a) argues
that the uncertainty inherent in dating individual glacial landforms is too high to make any
correlations with millennial-scale climate drivers. This perspective is further supported by Roe
and O’Neal (2009) who provide a case study proving that the stochastic variability of glacier
63

response to climate forcing can be misattributed to changes in climate. However, we can still
attain a synoptic understanding of paleoclimate dynamics using the current global dataset of
cosmogenic-nuclide exposure dated alpine glacier landforms (Balco, 2020a).
Correlation remains poorly quantified, i.e., lack of standardized statistical measures, in
the field of glacier chronology reconstruction via moraine records, most often only qualitatively
estimated through visual comparison of chronostratigraphic columns. Two situations can result
in false correlations between glacier response and climate events (Kirkbride and Winkler, 2012):
1. Type I (“lumping” errors): Multiple climatic events are correlated as a single event because
the glacier response time is less than resolution of the utilized dating methodology.
2. Type II (“splitting” errors): A single climatic event is distinguished into multiple events
because the timescales of glacier response do not match the chronological resolution, usually
caused by high-resolution dating methods.
Inadequately resolved glacier chronologies (i.e., greater dating uncertainty/poor statistical
confidence) have an increased likelihood of being "lumped" into regional-to-global correlations
because greater dating uncertainty allows moraines, which may have been generated from
multiple distinct climatic events, to appear as if they were generated from the same event
(Kirkbride and Winkler, 2012). Conversely, highly resolved glacier chronologies are more likely
to be "split" into many stages of glacier advance, making correlation difficult at broad spatial
scales. A narrow range of dating uncertainty may cause a singular climatic event to be
erroneously interpreted as numerous, distinguishable events (Kirkbride and Winkler, 2012).
Given the sparse availability of TCN dated surface exposure ages in the northeastern
Canadian Cordillera, I sourced 10Be ages from across WNA to make assumptions about the

64

climatic events influencing the Itsi Range. Consequently, the regional 10Be chronology likely
represents Type I “lumping” errors, as the gaussian distributions visualized in Figure 3.2
illustrate millennial-scale levels of 2σ uncertainty for each end group, far greater than the
average response time for alpine glaciers. This approach risks over-generalizing multiple
climatic events affecting glacial response, resulting in moraine formation, as representing a
singular climatic event. The benefit of Type I error lumping is that it simplifies the interpretation
of glacial records by reducing the number of discrete events to consider, potentially making it
easier to identify broad climatic trends and patterns over time. However, this simplification
comes at the cost of potentially overlooking the complexity and variability of past climate
events, which could lead to a misleading portrayal of glacial responses and climatic history.
3.5 Conclusion
I compare the glacier chronology at Itsi Glacier with 10Be moraine ages across WNA,
incorporating documented intervals of glacier advances compiled from multiple reviews using
radiocarbon dating, dendrochronology, and glaciolacustrine sedimentation records. To facilitate
comparisons, I categorize WNA 10Be moraine ages into four groups based on exposure age and
morphostratigraphic position relative to the glacier terminus.
The inner moraine at Itsi Glacier, dated to 0.20 ka ± 0.05 (n = 3), aligns well with the
regional glacier advance period at 0.40 – 0.06 ka. The compiled regional 10Be moraine record
yields an end group 1 age of 0.13 ka ± 0.10 with a MAD value of 67 years (Table 3.2). When
accounting for differences in latitude, northern latitudes show end group 1 moraines dating to
0.30 ka ± 0.10 (n = 11), while southern latitudes (50.9°) record emplacement at 0.09 ka ± 0.03 (n
= 14). This suggests that end group 1 moraine emplacement occurred decades to centuries earlier
at northern latitudes compared to southern latitudes.
65

The outer moraine at Itsi Glacier, dated to 0.96 ka ± 0.36 (n = 4), does not align with any
singular regional Neoglacial advance. Individually dated boulders accord to multiple advances:
2.09 ka ± 0.15 (228B04) at 2.20 – 1.80 ka advance, 1.16 ka ± 0.08 (228B01) at 1.71 – 1.20 ka
advance, and 0.43 ka ± 0.04 (228B05) at 0.90 – 0.40 ka advance. The regional 10Be moraine
record yields an end group 2 age of 0.47 ka ± 0.35 with a MAD value of 201 years (Table 3.2).
Although all end groups are non-parametric, end group 1 has a significantly lower MAD
value, indicating less geological scatter than end group 2. If accepting the proposed mechanism
of moraine amalgamation, Itsi Glacier sample 228B05 represents the true minimum-limiting age
for maximum LIA extents at 0.43 ka ± 0.04. According to the regional 10Be record and new ages
from the Itsi Glacier, Holocene maximum extents were reached centuries earlier in the
northeastern Canadian Cordillera when compared to alpine glaciers of the western and
southwestern Canadian Cordillera. Additional ages from erratic boulders suggest CIS melt
exposed alpine glaciers among the Itsi Range around 11.61 ka ± 0.54 (n = 2), consistent with
regional observations (Menounos et al., 2017).
Considering the effects of geologic scatter, particularly the roles of inheritance and postdepositional disturbance on Itsi Glacier moraines, the regional moraine record suggests that
inheritance on the outer moraine is more likely than ice-cored moraine induced post-depositional
disturbance. The most parsimonious explanation for consistent reports of 10Be geological scatter
on outer moraines from LIA moraine pairs in the northeastern Canadian Cordillera is that they
likely represent amalgamated moraines, heterogeneous mixtures of previously exposed and
reincorporated moraine boulders. Notably, the 10Be ages of end group 1 moraines throughout the
Selwyn and Mackenzie mountains exhibit significantly less scatter compared to the end group 2
moraines (Figure 3.2 and Figure 3.3). Moreover, the anomalous older boulder ages sampled on
66

end group 2 moraines align with pre-LIA moraines dated in AK and BC, centering around the
3.54 – 2.50 ka regional advance (Barclay et al., 2009; Menounos et al., 2013a; Pendleton et al.,
2017; Young et al., 2009).
Instead of obliterative overlap presiding, the outermost moraine might have posed a
topographical barrier high enough to prevent further glacier advance, leading to multiple
advances depositing boulders along the same ridge (Denton and Karlén, 1977; Pendleton et al.,
2017). Evidence for this mechanism, beyond intervals of regional advance correlative with
boulder ages, is lacking. Boulders exposed during less extensive glacial periods may simply be
recycled by the advancing glacier or inherited debris transported supraglacially before deposited
as a moraine (Balco, 2020a). More research quantifying the geomorphic uncertainties inherent to
TCND, e.g., Tomkins et al. (2021), must be done before declaring whether the geologic scatter is
recording “real” events or merely correlations between ages that happen to overlap with known
periods of advance.

4. Comparison of moraine record with climate proxy records
4.1 Introduction
Mountain glaciers act as global climate indicators due to their sensitivity to climate
forcings (Kaufman and Broadman, 2023). The presence of numerous periodicities (i.e., multiple
climate drivers at different frequencies) results in the reinforcement and modulation of climate
signals (Kirkbride and Winkler, 2012). The spatial variation of these climate signals may be
attributed to the complex interplay of atmospheric and oceanic circulation and teleconnection
patterns; therefore, globally synchronous climatic signals, such as a glacier advance, appear
when high-frequency cycles reinforce low-frequency cycles (Kirkbride and Winkler, 2012). The
67

Little Ice Age (LIA) represents an example of this, wherein minimal solar output (highfrequency; decadal-to-centennial scale), volcanic forcing (high-frequency; annual-to-decadal
scale), declining orbital insolation (low-frequency; 104-year timescale), and pre-industrial
decline of atmospheric CO2 combined to create a globally synchronous climatic signal
(Kirkbride and Winkler, 2012; Koch et al., 2019; Owens et al., 2017).
Quasi-periodic teleconnections, such as the Arctic Oscillation (AO), Pacific Decadal
Oscillation (PDO), and El Niño-Southern Oscillation (ENSO), demonstrably influence glacier
mass balance (Hodge et al., 1998; Moore et al., 2009); however, they are rarely detected in
moraine sequences due to the high-frequency (decadal) signal of these quasi-periodic oceanatmospheric circulation patterns (Kirkbride and Winkler, 2012). Because of this, attention is
given to low-frequency climate signals at the centennial-to-millennial scale.
4.2 Methodology
Fluctuations of glacier length can only distill an indirect, delayed, and filtered signal of
disturbances in mass balance due to the interplay of lag time components: reaction time and
response time (Kirkbride and Winkler, 2012; Oerlemans, 1994; Solomina et al., 2015).
Jóhannesson et al. (1989) define response time as “the time a glacier takes to adjust to all but a
factor of 1/e (37%) net change in mass balance,” making it a quantifiable constant for all glaciers
(Solomina et al., 2015). The equation for response time is:
𝑡𝑟 = 𝐻𝑀𝐴𝑋 /𝑎̇ 0

(14)

Where 𝑡𝑟 is response time (years), 𝐻𝑀𝐴𝑋 is mean thickness (m), and 𝑎̇ 0 is ablation rate near the
glacier terminus (m yr-1) (Cuffey and Paterson, 2010). Valley glacier response time generally
ranges from 10 – 50 years (Oerlemans, 1994). Reaction time, defined as the time lag between a
68

change in mass balance and response at the glacier’s terminus (Jóhannesson et al., 1989;
Solomina et al., 2015), allows for direct observation but cannot be quantified as a constant
among a suite of glaciers as it depends on a glacier’s history of mass balance change (Solomina
et al., 2015).
Response times vary significantly across alpine glaciers within the western North
America (WNA) 10Be chronology, but likely exist on the order of 101 – 102 years (Jóhannesson
et al., 1989). Using an estimated 27.3 m mean thickness of the Itsi Glacier (Farinotti et al., 2019),
and a mean yearly retreat rate of 1.93 m yr-1 from 1754 to 1949 CE estimated using the GTT tool
(Table Sx), 𝑡𝑟 is estimated to be 14 years. Given that the uncertainties for median ± interquartile
range (IQR) end group ages are on the centennial-to-millennial scale, the climate proxy
discussions are conducted at a scale where adjustments for 𝑡𝑟 is not required.
I qualitatively compare the developed glacier chronology to climate proxy records
including volcanic forcing, North Pacific Index (NPI), Total Solar Irradiance (TSI), GHG
emissions, orbital forcing, and temperature and precipitation anomalies. These datasets processed
using Python (v. 3.11.9) help determine if the proxy records support or refute assertions made by
terrestrial in situ cosmogenic nuclide (TCN) derived ages.
Volcanic forcing through Stratospheric Sulfur Injections (VSSI) significantly modulated
Earth's climate over the Holocene (Kaufman and Broadman, 2023; Sigl et al., 2022; Toohey and
Sigl, 2019). VSSI events of sufficient magnitude lead to a reduction of incoming solar insolation
at the Earth's surface, contributing to climate cooling at the decadal scale (Porter, 1986; Sigl et
al., 2022; Toohey and Sigl, 2019). A correlation coefficient of r = -0.6 between mean annual 2 m
air temperature and volcanic forcing in the Northern Hemisphere (NH) further supports this
relation (van Dijk et al., 2024). Increased clustering of volcanic eruptions can trigger ice-albedo
69

feedbacks and reduce ocean heat content, potentially inducing centennial scale cooling; the
cooling duration depends on the clustering, pacing, and magnitude of eruptions (Kaufman and
Broadman, 2023). Reconstructed volcanic eruption counts from Greenland and Antarctica icecore sulfate and sulfur records provide a globally resolved Holocene chronology of every
volcanic eruption with a mean uncertainty of ±1 – 5 years over the past 5 ka (Sigl et al., 2022).
Eruption counts only include those with a minimum of 5 Tg S injected into the stratosphere (Sigl
et al., 2021). The eVolv2k database provides supplementary VSSI values that record a 25%
reduction in VSSI over the 20th century as compared to a pre-industrial 2.45 – 0.05 ka mean
(Toohey and Sigl, 2019).
The NPI is associated with the Aleutian Low (AL) pressure centre which, a lowfrequency teleconnection pattern known to influence glaciation in British Columbia (BC) and
Alaska (AK) (Mood and Smith, 2015; Osterberg et al., 2014). Construction of the NPI proxy
involved measuring sodium concentrations by ion chromatography from an ice core collected at
the summit of Mt. Logan, YT (60.6° N, -140.5° W) assembled at a 3 year resolution (Osterberg
et al., 2014). Lower NPI values, indicating a stronger eastward AL shift, are associated with
increased precipitation and cooler temperatures conducive to glacial growth. Higher NPI values,
caused by westward AL shifts, bring dry and warm conditions, resulting in glacier retreat (Mood
and Smith, 2015; Osterberg et al., 2014).
Solar radiation emitted by the sun, TSI, controls the amount of solar forcing on Earth and
is reconstructed at a 5 year resolution by relating open solar magnetic field strength to 10Be
measurements in polar ice cores (Steinhilber et al., 2009; Wanner et al., 2008). TSI has been
demonstrated to influence mean surface temperature variations over the last millennium, with
prolonged negative TSI periods influencing changes in high-latitude NH glaciers (Gray et al.,
70

2010; Solomina et al., 2015). The most established mechanism explains that variations in solar
UV radiation directly affect the amount of ozone in the stratosphere, leading to temperature
changes, which then alter wind patterns and subsequently impact broader climate dynamics
(Gray et al., 2010). Amplification mechanisms, such as solar forcing triggering expansion of sea
ice, may also explain how minor TSI anomalies lead to glacier advances (Solomina et al., 2015).
However, there is no definite consensus concerning the relationship between solar variability and
Earth's climate system (Gray et al., 2010).
Köhler et al. (2017) compile annually resolved GHG data from 13 globally distributed ice
cores and atmospheric observations. While the CO2 proxy reflects this global compilation, CH4
data are collected solely from the SH to account for discrepancies in the timing of major climatic
events between hemispheres (Köhler et al., 2017). This CH4 data is estimated to represent the
lower boundary of global CH4 concentrations (Köhler et al., 2017). As GHG levels rise, they
positively contribute to radiative forcing, directly leading to increased global temperatures
(Köhler et al., 2017; Tierney et al., 2020). CO2 is the primary contributor, accounting for twothirds of radiative forcing,
𝐶𝑂2
)
Δ𝑅[𝐶𝑂2] = 5.35 ⋅ ln (
278 𝑝𝑝𝑚

(15)

while CH4 plays a smaller, but still significant, role (Köhler et al., 2017).
Δ𝑅[CH4] ∼ 1.4 ⋅ 0.036 ⋅ (√CH4 − √742 𝑝𝑝𝑏)

(16)

Both GHGs are converted to and plotted as radiative forcing (W m-2) using equation (15 for CO2
and equation (16 for CH4, where CO2 is measured concentration in parts per million (ppm) and
CH4 is measured concentration in parts per billion (ppb) (Köhler et al., 2017; Myhre et al., 1998).
71

Orbital variations, or Milankovitch cycles, significantly influence Earth’s glacialinterglacial cycles, earning the colloquialism “pacemaker of the ice ages” (Berger, 1988;
Marshall, 2011). These cycles, consisting of eccentricity, obliquity, and precession, affect solar
insolation distribution and align cyclically on 104 – 105-year scales, leading to cooler or warmer
summers at high latitudes (Marshall, 2011; Wanner et al., 2008). NH orbital forcing primarily
controls glacial-interglacial cycles by altering solar insolation distribution over the Earth (Berger,
1988; Marshall, 2011). A decrease in NH solar insolation at high latitudes can trigger glacial
periods due to ice-albedo feedback loops that result from increased cryospheric albedo (Kaufman
and Broadman, 2023). By 1950, summer insolation at 65° N decreased by 37.91 W m-2 relative
to the beginning of the Holocene at 12 ka and continues to steadily decline (Berger and Loutre,
1991).
I obtained temperature and precipitation anomalies from the TraCE-21ka transient
climate simulation (Liu et al., 2009), centered on the Itsi Range (61.2° N, -131.2° W) with a
temporal resolution of 10 years and a spatial resolution of 3.68° latitude x 3.75° longitude. Both
variables use 5 – 3 ka as the reference period for anomaly calculations.

72

4.3 Results and discussion

Figure 4.1: Neoglacial climate forcings and indicators. (a) Global volcanic eruptions per decade
(30-year moving average) and major VSSI event magnitudes (Tg S) (Sigl et al., 2021; Toohey
and Sigl, 2019). (b) NPI at Mt. Logan (750-year spline) (Osterberg et al., 2014). (c) TSI
73

anomalies (ΔW m-2) plotted as 50-year moving averages (light line) and smoothed using a
Gaussian filter (dark line) (Steinhilber et al., 2009). (d) A spline smoothed dataset of atmospheric
GHGs CO2 and CH4 plotted as radiative forcing (ΔW m-2) (Köhler et al., 2017) Increases in
GHG concentrations past 0.1 ka are not plotted. (e) July mean Solar Insolation at 65° N (W m-2)
at a 1000 year resolution (Berger and Loutre, 1991). (f) Itsi Range temperature (ΔT) and
precipitation anomalies (ΔP) (Liu et al., 2009) plotted as 50-year moving averages. (g) Regional
compilation of late Holocene 10Be boulder ages, including those collected for this thesis (Dortch
et al., 2010; Hawkins et al., 2023, 2021; Howley, 2008). Blue and orange lines represent the sum
of each 2σ Gaussian distribution for end groups 1 and 2, respectively. Itsi Glacier moraine end
groups are highlighted separately within the distributions.
Analyzing the terrestrial record of 10Be ages to regionally relevant climate forcings can
help deduce which climate forcings influence glacier expansion and retreat in WNA. Rather than
evaluating each proxy individually, I focus on broader patterns at multicentennial-to-millennial
scales by analyzing moraine end groups. This approach accounts for the high spatiotemporal
variability among proxy types and acknowledges that multiple proxies may collectively
influence moraine deposition.
Decreasing solar insolation primarily drove the steadily declining temperature anomalies
observed in the Itsi Range, responsible for LIA glacial expansion to Holocene maximum extents
in the NH (Clague et al., 2009; Kaufman and Broadman, 2023; Menounos et al., 2009; Wanner
et al., 2011). While the onset of Neoglacial glacier growth varies widely across WNA, highfrequency climate forcings can obscure long-term low-frequency trends but still provide valuable
insights into regional climate changes (Solomina et al., 2015). Understanding these forcings
helps elucidate the mechanisms driving moraine deposition and glacier advances in WNA
throughout the late Holocene.
4.3.1 Moraine end groups
Median end group 2 regional ages at 0.47 ka ± 0.35 (n = 36) variably correspond with the
climate proxies, while local Itsi Glacier ages at 0.96 ka ± 0.36 (n = 4) have weaker correlations

74

(Figure 4.1). A VSSI spike of 5.33 Tg S at 2.4 ka coincides with a negative TSI anomaly,
preceding a period of cooling that began at 2.3 ka and culminated in the 2.2 – 1.80 ka regional
advance. Sample 228B04, dated at 2.10 ka ± 0.15, may reflect this period of cooling if its
deposition age is reliable.
While GHG levels remained relatively stable throughout the mid-Holocene, a sustained
increase occurred from 1.0 – 0.4 ka, possibly linked to the agricultural practices of pre-industrial
civilizations (Ruddiman et al., 2016). This rise led to an increase in radiative forcing by about
0.46 W m-2 compared to the mean at 6 ka (Kaufman and Broadman, 2023), possibly
counteracting natural cooling expected from orbital cycles. However, the impact would likely
have been subtle compared to the more pronounced effects of natural climate variability (Roe
and O’Neal, 2009).
A TSI minimum influencing glacial advance around 0.5 ka precedes the ages of many
regional end group 2 moraines. Although the NPI is loosely reflected in temperature trends
overall, an increase from 0.6 – 0.4 ka and the subsequent temperature rise from -0.3°C to -0.1°C
appears to have facilitated end group 2 moraine deposition. However, the local Itsi Glacier end
group 2 shows limited correlation with all climate proxies.
Median end group 1 regional ages at 0.13 ka ± 0.10 (n = 25) and local Itsi Glacier ages at
0.20 ka ± 0.05 (n = 3) show concordance with all assessed climate proxies (Figure 4.1). Volcanic
aerosol injections have likely induced global cooling over the last two millennia (Kaufman and
Broadman, 2023), with volcanic forcing contributing significantly to the high natural climate
variability throughout the LIA (van Dijk et al., 2024). Increased clustering of eruptions starting
around 1.3 ka likely amplified orbitally forced Neoglacial cooling, culminating in the LIA, as
supported by Holocene transient simulations (Bader et al., 2020; van Dijk et al., 2024).
75

TSI and NPI trends correlate well with end group 1 moraine deposition. TSI anomalies
below zero at 0.7 ka coincide with the onset of the LIA, persisting until 0.2 ka when positive
anomalies emerged. This corresponds with the culmination of regional end group 1 ages.
Decreasing NPI values from 0.4 ka align with increasing precipitation and decreasing
temperature trends, facilitating the deposition of most end group 1 moraines.
From 0.4 – 0.3 ka, global CO2 emissions gradually dropped by 0.5 ppm per decade
following significant land-use changes in the Americas following contact between the New
World and the Old World (King et al., 2024; Koch et al., 2019). Anthropogenic GHG emissions,
starting around 0.2 ka, may have prematurely ended Neoglaciation as CO2 and CH4 levels began
to rise (Wanner et al., 2011).
4.3.2 Anthropogenic warming
Although many NH glaciers were nearly as small or smaller than current extents from 9 –
5.5 ka, the observed global retreat of glaciers in the last century does not align with expected
orbital forcing patterns expected to encourage glaciation and is primarily attributed to
anthropogenic influences (Marzeion et al., 2014; Solomina et al., 2015). Agassiz Ice Cap core
data reveal that current melt rates are the highest in 4200 years, comparable to the early
Holocene thermal maximum (Fisher et al., 2012). Between 1851 and 2010 CE, approximately 25
± 35% of global glacier mass loss is attributed to anthropogenic warming, a figure that rises to 69
± 24% from 1991 to 2010 in western Canada (Marzeion et al., 2014).
Arctic amplification, driven by NH sea ice expansion, results in Arctic summer
temperature variations 3 – 4 times larger than global changes, explaining the sensitivity of early
Holocene and modern Arctic glaciers (Miller et al., 2013, 2010). Century-scale early Holocene

76

advances thought to be caused by cooling from LIS meltwater, which weakened the Atlantic
thermohaline circulation, mirror modern concerns with Greenland ice sheet melt (Menounos et
al., 2017; Nesje et al., 2004).

Figure 4.2: Projected Itsi Range (62.9°N, -130.2°W,) SSP climate change scenarios (Cannon et
al., 2015; ClimateData.ca, 2024; McKenney et al., 2011; Tam et al., 2019).
Shared Socioeconomic Pathways (SSPs) paired with different levels of radiative forcing
(W m-2), measured as positive deviations from a pre-industrial baseline, assess GHG influence
on the energy balance of the climate system (Riahi et al., 2017) (Figure 4.2). COP26 climate
pledges indicate an estimated global mean temperature increase of +2.7°C, corresponding to
SSP2-4.5, which would result in near or total deglaciation of WNA glaciers (95 ± 8%) by the
2100 CE (Rounce et al., 2023).
4.4 Limitations
While correlating glacier chronologies with climate proxies can validate paleoclimatic
patterns, it does not imply causation. The absence of formal quantitative criteria for correlating
77

local glacier chronologies means that only qualitative conclusions can be drawn (Kirkbride and
Winkler, 2012). Consequently, the interpretation of these proxy records in relation to moraine
deposition is inherently subjective.
The lack of quantification for the response times of the glaciers represented by each end
group introduces additional uncertainty when interpreting their relationship to climate proxies.
The assumptions made about the moraines belonging to different end groups directly affect the
interpretation of their Gaussian distributions in relation to the proxy records. Different
approaches to grouping or analyzing the moraines could lead to alternative interpretations of the
same data.
Proxies each have their own inherent limitations. Volcanic eruptions and associated VSSI
aerosol spread have high spatiotemporal variability, meaning that it cannot be stated that all
eruptions affected NH climate at high latitudes in the northeastern Canadian Cordillera. The
NPI’s 750-year low-resolution spline simplifies climatic variations at fine temporal scales,
offering only a broad approximation helpful for relating low-frequency proxies (Osterberg et al.,
2014). The correlation between 2 m air temperature and TSI is weak (r = 0.06), except during
major LIA solar minima events at 0.56 – 0.39 ka and 0.33 – 0.23 ka. (van Dijk et al., 2024).
The Mt. Logan ice core record as reported by Osterberg et al. (2014) may be subject to
revision. Preliminary results from a 2022 Mt. Logan ice core study (Holland et al., 2023) suggest
significantly different age estimates (Holland et al., 2023). However, no formal publication of
these findings has been released to date.
Other factors on glaciation in the northeastern Canadian Cordillera were not addressed
due to limited proxy records or the complications associated with correlating high-frequency

78

forcings at multicentennial-to-millennial scales. These include internal climate system variability
(Roe and O’Neal, 2009) and the PDO known to influence western Canadian glacier behavior
(Koch et al., 2009). When considering all these limitations, the qualitative correlations drawn
between the glacier chronology and these proxies should be interpreted with caution.
4.5 Conclusion
A reduction of orbital insolation and increased volcanic forcing initiated LIA glaciation
(Miller et al., 2012; van Dijk et al., 2024). The regional compilation of late Holocene 10Be ages
shows two distinct peaks in glacial advances during this period, with median ± IQR ages of
WNA moraine deposition at 0.47 ka ± 0.35 and 0.13 ka ± 0.10. These peaks show some
alignment with the climate proxies in Figure 4.1, indicating a synchronous regional response to
certain multicentennial-scale climate forcings, many other peaks lack a coherent counterpart.
Both regional end group ages and Itsi Glacier end group 1 show stronger coherence with
the climate proxies, in contrast to Itsi Glacier end group 2 ages which does not correlate well
with any proxy beyond orbital forcing. This lack of synchronicity likely indicates significant
scatter in the Itsi Glacier end group 2 ages, deviating from the regional 10Be age signature.
Apart from orbital forcing, proxies that most closely align with the 10Be dated end group
1 ages – and likely exert the most influence on glacial advance and retreat – include temperature
and precipitation anomalies, GHG emissions, and TSI anomalies. Proxies with a less direct
impact on end group 1 moraine deposition are NPI and volcanic forcing. For regional end group
2, TSI anomalies show the strongest correspondence, although temperature anomalies, GHG
emissions, and NPI display plausible but less direct correlations when accounting for
uncertainties.

79

The Mt. Logan NPI values suggest shifts in the AL significantly influenced glacial
events. A westward shift of the AL likely drove glacial retreat after the deposition of end group
2. A subsequent eastward shift after end group 2's deposition contributed to another, though less
extensive, glacial advance that culminated in end group 1 deposition. The rise in GHG emissions
appears to have prematurely ended the advances of end group 1, despite NPI values, orbital
trends, and temperature and precipitation data indicating conditions favorable for glacier growth.
This suggests that anthropogenic influences have increasingly disrupted the natural drivers of
glacier advance and retreat in this region.
The timing of end group 1 and 2 moraine deposition aligns with periods of negative TSI
anomalies, stronger AL influence, and reduced summer insolation, supporting the possibility that
these proxies collectively drove moraine deposition in the northeastern Canadian Cordillera.

5. Conclusions
5.1 Summary of advances
This thesis presents a latest Pleistocene and Holocene chronology for the Itsi Glacier in
the Selwyn Mountains, YT using 10Be surface exposure dating of moraine boulders to determine
onset of glacial retreat. This research marks the first application of cosmogenic nuclide dating on
alpine cirque glacier moraines within the Itsi Range, YT. Erratics immediately beyond the Itsi
Glacier terminal moraine date to 11.61 ka ± 0.54 (n = 2) and indicate that the CIS was no more
extensive than the erratics’ position at this time. The Itsi Glacier reached its maximum Holocene
extent at 0.96 ka ± 0.36 (n = 4) as recorded by the outer moraine. Considering the proposed
method of moraine amalgamation, I interpret Itsi Glacier sample 228B05 as the true minimumlimiting age at 0.43 ka ± 0.04 as related to the regional advance between 0.90 – 0.40 ka and
80

aligning with the timing of other maximum Holocene extents in the Selwyn and Mackenzie
Mountains as concluded by Hawkins et al. (2023). These findings suggest that LIA maximum
extents were reached earlier (~0.40 ka) in the northeastern Canadian Cordillera than the 19thcentury maxima recorded in the western and southwestern Cordillera (Menounos et al., 2009;
Mood and Smith, 2015). Following the Holocene maximum, the Itsi Glacier experienced a minor
re-advance at 0.20 ka ± 0.05 (n = 3), as recorded by the inner moraine. This date aligns with the
regional compilation of Neoglacial advance intervals in WNA from 0.40 – 0.06 ka discussed in
Chapter 3.3.2.
Comparing end group 1 and end group 2 moraines in the Selwyn and Mackenzie
Mountains reveals patterns of geologic scatter. All end group 1 moraines exhibit significantly
less age scatter, in terms of NMAD and X2, and stronger correlations with regional climate
proxies, suggesting a coherent glacial signal. Latest Holocene advances are primarily driven by
orbital forcing, with other key influences including temperature and precipitation anomalies,
GHG emissions, and TSI fluctuations.
In contrast, all end group 2 moraines show low concordance with climate proxies. This
lack of synchronicity likely reflects significant geologic scatter, deviating from the expected 10Be
age signature for the region. Anomalously old 10Be ages from end group 2 moraines across the
Selwyn and Mackenzie Mountains (2.29 ka ± 0.83, n = 10) appear correlative with end group 3
moraines at northern latitudes in Alaska (2.33 ka ± 0.63, n = 20) (Pendleton et al., 2017; Young
et al., 2009) and an end group 3 moraine at Tiedemann Glacier in British Columbia (2.60 ± 0.28,
n = 6) (Menounos et al., 2013a). These end group 3 moraines are coeval to the Tiedemann-Peyto
Advance between 3.54 – 2.50 ka. The consistent geologic scatter in end group 2 likely reflects

81

amalgamated moraines, heterogeneous mixtures of previously exposed and reincorporated
moraine boulders, deposited during the Tiedemann-Peyto Advance.
Photogrammetric analysis of the Itsi Range glaciers reveals an accelerating mass loss.
Between 1949 and 2004, glaciers in the region experienced a geodetic mass balance loss of
-0.190 ± 0.065 m w.e. yr-1. This rate more than doubled to -0.382 ± 0.136 m w.e. yr-1 during the
period 2004 to 2022, representing a 101% increase in mass loss.
5.2 Research limitations
A key challenge when conducting TCND of alpine glaciers is managing the
methodological issue of geologic scatter. While the Itsi Glacier’s inner moraine and outboard
erratics exhibit minimal scatter of 10Be ages, the outer moraine typifies this recurring problem.
Of the four boulders sampled from the outer moraine, two yielded ages from the past
millennium, while the other two indicated 1 – 2 ka of surface exposure. Presence of geologic
scatter complicates the attribution of a moraine to a specific glacial advance, as the true age is
obscured by these outliers.
The high cost of remote sampling and laboratory processing of boulders for 10Be
concentration limits the ability to thoroughly assess late Holocene fluctuations at the Itsi Glacier.
Collecting three to five samples per moraine is common given the expense, but this sample size
often proves insufficient when outliers are present, as seen with the Itsi Glacier’s outer moraine,
reducing the reliability of the dataset for robust conclusions.
Using median ± IQR is advantageous because it does not assume normality of data and is
less sensitive to outliers, such as those found in moraines with scattered boulder ages. However,
this approach can result in less certainty about dispersion/variance compared to measures like
82

mean ± standard deviation, leading to broader error margins. These larger error bounds reduce
the chronological resolution of moraine ages, increasing the risk of overinterpretation due to the
increased apparent similarity of moraine ages; a challenge specifically when attributing the Itsi
Glacier outer moraine to a regional advance as discussed in Chapter 3.
The geomorphic uncertainties inherent to TCND limit the ability to confidently attribute
geologic scatter to moraine amalgamation. While the alignment of anomalously old boulder ages
with pre-LIA advances suggests amalgamation as a possible mechanism, other factors, namely
inheritance or post-depositional disturbance, cannot be ruled out. Without further research
quantifying these uncertainties, it remains difficult to assert whether moraine amalgamation is
the definitive cause of the scatter.
The Chapter 4 discussion of moraine ages in relation to climate forcings is also
vulnerable to overinterpretation. This is because it relies on subjective visual correlations among
climate proxies, many of which have large uncertainties. As a result, interpretation of climate
forcings as related to moraine ages may only suggest broad, correlated trends, rather than
definitive evidence for causation of specific glacial advances.
5.3 Research recommendations
1. Utilize tandem dating methods for geochronological research. Both 14C and 10Be can be
extracted from quartz, enabling paired 14C-10Be analysis a way to build a more robust
chronology (Lamp et al., 2019). This technique leverages differences in half-life decay
and nuclide production to detect disequilibrium in exposure ages from the LGM to the
late Holocene (Hippe, 2017; Schimmelpfennig et al., 2022). By identifying inheritance or
exhumation processes that contribute to geologic scatter on moraines, paired dating

83

improves the reliability of exposure ages (Goehring et al., 2022; Stroeven et al., 2014). A
more cost-effective option is utilizing organic material when available. Radiocarbon
dating from glacier forefields can complement 10Be dating, improving the precision of
glacial chronologies and constraining the duration of Holocene advances at individual
sites (Hawkins et al., 2021). However, obtaining organic 14C in glacier forefields is often
not possible because many sites at far northern latitudes, including the glaciers of the Itsi
Range, lack the growing seasons necessary to support organic growth such as trees.
2. Strengthen understanding of geomorphic context in relation to moraine boulder
exposure ages. While sampling criteria were carefully considered in this study,
uncertainties remain regarding best practices for sampling moraines as indicated by
emerging research from Tomkins et al. (2021). For latest Holocene research that targets
younger moraines (<1 ka), boulder sampling should avoid moraine crests and instead
target more stable intra-landform locations on ice-proximal and ice-distal slopes
(Tomkins et al., 2021). These preliminary conclusions about sampling strategies need
further quantitative testing, with the framework of tandem TCND with Schmidt hammer
sampling proposed by Tomkins et al. (2021, 2018) a viable path forward.
3. Further research mechanisms for moraine amalgamation. Although uncertainties for
measuring concentrations 10Be have steadily decreased over the past thirty years,
observations of geologic scatter remain reflecting geomorphic uncertainties in the field
(Schaefer et al., 2022). Developing methods to distinguish amalgamated moraines from
those representing single glacier advances is crucial for resolving persistent
methodological issues of geologic scatter. One approach could compare X2 or NMAD
values across a large dataset of 10Be-dated moraines grouped by morphostratigraphic

84

position. Another strategy involves resampling previously dated moraines with signs of
geologic scatter and correlating these to well-constrained local advances to qualitatively
assess processes of moraine amalgamation in a more controlled manner.

Appendix

Figure A.1: Regional compilation of WNA late Holocene 10Be boulder ages, including those
collected for this thesis, sampled from moraines (Dortch et al., 2010; Hawkins et al., 2023, 2021;
Howley, 2008; Menounos et al., 2013a; Pendleton et al., 2017; Young et al., 2009). Thin grey
curves show ages of individual boulders rendered as Gaussian distributions. The thick curve
represents the sum of each Gaussian distribution for the entire population of 10Be moraine end
groups 1 – 3 ages. Itsi Glacier moraine end groups 1 and 2 are inclusive in the orange curve and
also highlighted separately within the summed Gaussian distributions. Corresponding box-andwhisker plots denote the median ± IQR, minimum, and maximum surface exposure ages of each
moraine group.

85

Figure A.2: Regional compilation of WNA late Holocene scatter plots showing 10Be moraine
ages, including those collected for this thesis, grouped by morphostratigraphic position (Dortch
et al., 2010; Hawkins et al., 2023, 2021; Howley, 2008; Menounos et al., 2013a; Pendleton et al.,
2017; Young et al., 2009). Y-axis displays IQR (ka) and x-axis displays n, latitude, and
longitude.

86

Figure A.3: Output from the GTT toolbox showing line connections between traced Itsi Range
glacier termini and corresponding moraines used in retreat rate calculations.

Figure A.4: Compilation of all sampled boulders with their associated sample number. Photos
taken by Gerald Osborn.

87

Figure A.5: Compilation of Itsi Glacier moraine field photos. Red arrows indicate true north.
Photos taken by Brian Menounos.

88

Figure A.6: Maps of mean elevation change among Itsi Range glaciers between the periods 1949
– 2004 and 2004 – 2022 for geodetic mass balance calculations.
Table A.1: Intervals of major WNA advances discussed in chapter 3 and illustrated in figure 3.1.
89

Advance
Locationsa
Evidenceb
References
(ka BP)
Early-middle Neoglacial glacier advances (4.40 – 3.40 ka)
4.40 – 3.97 CM, WCC
RC
(Menounos et al., 2009; Mood and Smith, 2015)
3.70 – 3.40 CM
RC
(Mood and Smith, 2015)
Middle-late Neoglacial glacier advances (3.54 – 1.20 ka)
3.54 – 2.50 AK, CM,
D, RC
(Barclay et al., 2009; Menounos et al., 2009;
WCC
Mood and Smith, 2015)
10
2.20 – 1.80 AK, CM
D, RC, Be (Barclay et al., 2009; Hawkins et al., 2021)
1.71 – 1.20 AK, CM,
D, RC, 10Be (Barclay et al., 2009; Menounos et al., 2009;
WCC
Mood and Smith, 2015)
Latest Neoglacial glacier advances (0.90 – 0.06 ka)
0.90 – 0.40 AK, CM, SM, D, RC, GL, (Barclay et al., 2009; Hawkins et al., 2021;
10
SN, WCC
Be
Menounos et al., 2009; Mood and Smith, 2015;
Solomina et al., 2015; Tomkins et al., 2008)
0.50 – 0.24 AK, SM,
D, RC, GL (Barclay et al., 2009; Menounos et al., 2009;
WCC
Tomkins et al., 2008)
0.40 – 0.06 AK, CM, SM, D, RC, GL, (Barclay et al., 2009; Hawkins et al., 2021;
10
SN, WCC
Be
Menounos et al., 2009; Mood and Smith, 2015;
Solomina et al., 2015; Tomkins et al., 2008)
a

Locations:
AK = Alaska, USA
CM = Coast Mountains, British Columbia, Canada
SM = Selwyn Mountains, Yukon Territory, Canada
WCC = Western Canadian Cordillera (encompassing Yukon Territory, Northwest Territories, British Columbia, and
Alberta)
SN = Sierra Nevada Mountains, USA
b
Evidence:
D = Dendrochronology
RC = Radiocarbon dating
10
Be = 10Be terrestrial in situ cosmogenic nuclide dating
GL = Glaciolacustrine sedimentation records

Table A.2: Change in area among all Itsi Range glaciers observed from historical orthophotos
(1949 and 2002) and satellite imagery (2022).
Total area (km2)
Year

Glacier
count

Max

Median

Min

Standard
deviation

Average

Sum

ΔArea
%

1949

18

2.89

0.97

0.09

0.78

1.08

19.51

-

2004

22

2.56

0.42

0.06

0.66

0.66

14.44

-26%

2022

25

1.99

0.22

0.03

0.53

0.45

11.17

-23%

90

References
Bader, J., Jungclaus, J., Krivova, N., Lorenz, S., Maycock, A., Raddatz, T., Schmidt, H., Toohey,
M., Wu, C.-J., Claussen, M., 2020. Global temperature modes shed light on the Holocene
temperature conundrum. Nat. Commun. 11, 4726.
Balco, G., 2020a. Glacier change and paleoclimate applications of cosmogenic-nuclide exposure
dating. Annu. Rev. Earth Planet. Sci. 48, 21–48.
Balco, G., 2020b. A prototype transparent-middle-layer data management and analysis
infrastructure for cosmogenic-nuclide exposure dating. Geochronology 2, 169–175.
Balco, G., 2011. Contributions and unrealized potential contributions of cosmogenic-nuclide
exposure dating to glacier chronology, 1990–2010. Quat. Sci. Rev. 30, 3–27.
Barclay, D.J., Wiles, G.C., Calkin, P.E., 2009. Holocene glacier fluctuations in Alaska. Quat.
Sci. Rev. 28, 2034–2048.
Barr, I.D., Lovell, H., 2014. A review of topographic controls on moraine distribution.
Geomorphology 226, 44–64.
Beedle, M.J., Menounos, B., Wheate, R., 2014. An evaluation of mass-balance methods applied
to Castle Creek Glacier, British Columbia, Canada. J. Glaciol. 60, 262–276.
Berger, A., 1988. Milankovitch theory and climate. Rev. Geophys. 26, 624–657.
Berger, A., Loutre, M.-F., 1991. Insolation values for the climate of the last 10 million years.
Quat. Sci. Rev. 10, 297–317.
Berger, A.L., 1978. Long-Term Variations of Caloric Insolation Resulting from the Earth’s
Orbital Elements1. Quat. Res. 9, 139–167.
Bierman, P.R., Bender, A.M., Christ, A.J., Corbett, L.B., Halsted, C.T., Portenga, E.W.,
Schmidt, A.H., 2021. Dating by cosmogenic nuclides, in: Encyclopedia of Geology.
Academic Press Oxford, pp. 101–115.
Bonsal, B., Shrestha, R.R., Dibike, Y., Peters, D.L., Spence, C., Mudryk, L., Yang, D., 2020.
Western Canadian freshwater availability: Current and future vulnerabilities. Environ.
Rev. 28, 528–545.
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.
Braumann, S.M., Schaefer, J.M., Neuhuber, S.M., Lüthgens, C., Hidy, A.J., Fiebig, M., 2021.
Early Holocene cold snaps and their expression in the moraine record of the eastern
European Alps. Clim. Past 17, 2451–2479.
Braumann, S.M., Schaefer, J.M., Neuhuber, S.M., Reitner, J.M., Lüthgens, C., Fiebig, M., 2020.
Holocene glacier change in the Silvretta Massif (Austrian Alps) constrained by a new
10Be chronology, historical records and modern observations. Quat. Sci. Rev. 245.
Broecker, W.S., Denton, G.H., Edwards, R.L., Cheng, H., Alley, R.B., Putnam, A.E., 2010.
Putting the Younger Dryas cold event into context. Quat. Sci. Rev. 29, 1078–1081.
Brown, L., Klein, J., Middleton, R., Sacks, I.S., Tera, F., 1982. 10 Be in island-arc volcanoes and
implications for subduction. Nature 299, 718–720.
Cannon, A.J., Sobie, S.R., Murdock, T.Q., 2015. Bias correction of GCM precipitation by
quantile mapping: how well do methods preserve changes in quantiles and extremes? J.
Clim. 28, 6938–6959.

91

Carlson, A.E., LeGrande, A.N., Oppo, D.W., Came, R.E., Schmidt, G.A., Anslow, F.S.,
Licciardi, J.M., Obbink, E.A., 2008. Rapid early Holocene deglaciation of the Laurentide
ice sheet. Nat. Geosci. 1, 620–624.
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.
ClimateData.ca, 2024.
Cogley, J., Hock, R., Rasmussen, L., Arendt, A., Bauder, A., Braithwaite, R., Jansson, P., Kaser,
G., Möller, M., Nicholson, L., Zemp, M., 2011. Glossary of glacier mass balance and
related terms. https://doi.org/10.5167/uzh-53475
Colpron, M., 2022. Yukon Bedrock Geology Map.
contributors, xdem, 2024. xdem. https://doi.org/10.5281/zenodo.11204531
Cook, E.R., Kairiukstis, L.A., 1990. Methods of dendrochronology: applications in the
environmental sciences. Kluwer Academic Publishers.
Crump, S.E., Anderson, L.S., Miller, G.H., Anderson, R.S., 2017. Interpreting exposure ages
from ice-cored moraines: a Neoglacial case study on Baffin Island, Arctic Canada. J.
Quat. Sci. 32, 1049–1062.
Cuffey, K.M., Paterson, W.S.B., 2010. The physics of glaciers. Academic Press.
Dalton, A.S., Dulfer, H.E., Margold, M., Heyman, J., Clague, J.J., Froese, D.G., Gauthier, M.S.,
Hughes, A.L., Jennings, C.E., Norris, S.L., 2023. Deglaciation of the north American ice
sheet complex in calendar years based on a comprehensive database of chronological
data: NADI-1. Quat. Sci. Rev. 321, 108345.
Darvill, C., Menounos, B., Goehring, B., Lesnek, A., 2022. Cordilleran Ice Sheet Stability
During the Last Deglaciation. Geophys. Res. Lett. 49, e2021GL097191.
Darvill, C., Menounos, B., Goehring, B., Lian, O.B., Caffee, M., 2018. Retreat of the western
Cordilleran Ice Sheet margin during the last deglaciation. Geophys. Res. Lett. 45, 9710–
9720.
Darvill, C.M., 2013. Cosmogenic nuclide analysis. Geomorphol. Tech. Br. Soc. Geomorphol.
Lond. 1–25.
Davis Jr, R., Schaeffer, O.A., 1955. Chlorine‐36 in nature. Ann. N. Y. Acad. Sci. 62, 107–121.
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.
Dommenget, D., 2009. The ocean’s role in continental climate variability and change. J. Clim.
22, 4939–4952.
Dortch, J.M., Owen, L.A., Caffee, M.W., Brease, P., 2010. Late Quaternary glaciation and
equilibrium line altitude variations of the McKinley River region, central Alaska Range.
Boreas 39, 233–246.
Erb, M.P., McKay, N.P., Steiger, N., Dee, S., Hancock, C., Ivanovic, R.F., Gregoire, L.J.,
Valdes, P., 2022. Reconstructing Holocene temperatures in time and space using
paleoclimate data assimilation. Clim Past 18, 2599–2629. https://doi.org/10.5194/cp-182599-2022
Evans, D.J., 2009. Controlled moraines: origins, characteristics and palaeoglaciological
implications. Quat. Sci. Rev. 28, 183–208.
Farinotti, D., Huss, M., Fürst, J.J., Landmann, J., Machguth, H., Maussion, F., Pandit, A., 2019.
A consensus estimate for the ice thickness distribution of all glaciers on Earth. Nat.
Geosci. 12, 168–173.

92

Fick, S.E., Hijmans, R.J., 2017. WorldClim 2: new 1‐km spatial resolution climate surfaces for
global land areas. Int. J. Climatol. 37, 4302–4315.
Fisher, D., Zheng, J., Burgess, D., Zdanowicz, C., Kinnard, C., Sharp, M., Bourgeois, J., 2012.
Recent melt rates of Canadian arctic ice caps are the highest in four millennia. Glob.
Planet. Change 84, 3–7.
Gibbons, A.B., Megeath, J.D., Pierce, K.L., 1984. Probability of moraine survival in a succession
of glacial advances. Geology 12, 327–330.
Girod, L., Nuth, C., Kääb, A., McNabb, R., Galland, O., 2017. MMASTER: improved ASTER
DEMs for elevation change monitoring. Remote Sens. 9, 704.
Goehring, B., Menounos, B., Osbron, 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.
Gordey, S.P., 2013. Evolution of the Selwyn Basin region, Sheldon Lake and Tay River map
areas, central Yukon. Geol. Surv. Can. Bulletin 599, 176. https://doi.org/10.4095/293034
Gosse, J.C., Phillips, F.M., 2001. Terrestrial in situ cosmogenic nuclides: theory and application.
Quat. Sci. Rev. 20, 1475–1560.
Gray, L.J., Beer, J., Geller, M., Haigh, J.D., Lockwood, M., Matthes, K., Cubasch, U.,
Fleitmann, D., Harrison, G., Hood, L., 2010. Solar influences on climate. Rev. Geophys.
48.
Groisman, P.Y., Karl, T.R., Knight, R.W., 1994. Observed impact of snow cover on the heat
balance and the rise of continental spring temperatures. Science 263, 198–200.
Hawkins, A.C., 2023. Multiple geochronological approaches to constrain late Quaternary glacial
histories in western Canada.
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.
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.
Heisinger, B., Lal, D., Jull, A., Kubik, P., Ivy-Ochs, S., Neumaier, S., Knie, K., Lazarev, V.,
Nolte, E., 2002. Production of selected cosmogenic radionuclides by muons: 1. Fast
muons. Earth Planet. Sci. Lett. 200, 345–355.
Heyman, J., Stroeven, A.P., Harbor, J.M., Caffee, M.W., 2011. Too young or too old: evaluating
cosmogenic exposure dating based on an analysis of compiled boulder exposure ages.
Earth Planet. Sci. Lett. 302, 71–80.
Hippe, K., 2017. Constraining processes of landscape change with combined in situ cosmogenic
14C-10Be analysis. Quat. Sci. Rev. 173, 1–19.
Hodge, S.M., Trabant, D.C., Krimmel, R.M., Heinrichs, T.A., March, R.S., Josberger, E.G.,
1998. Climate variations and changes in mass of three glaciers in western North America.
J. Clim. 11, 2161–2179.
Holland, K.M., Criscitiello, A.S., McConnell, J., Winski, D., Markle, B., Campbell, S.W.,
Kreutz, K.J., Chellman, N., Wensman, S., McCrimmon, D., 2023. Age scale and
preliminary findings of the 2022 Mt. Logan ice core, southwest Yukon. Presented at the
AGU Fall Meeting Abstracts, pp. C43A-07.

93

Howley, M.W., 2008. A late glacial and Holocene chronology of the Castner Glacier, Delta
River valley, Alaska. University of New Hampshire.
Hugonnet, R., Brun, F., Berthier, E., Dehecq, A., Mannerfelt, E.S., Eckert, N., Farinotti, D.,
2022. Uncertainty analysis of digital elevation models by spatial inference from stable
terrain. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 15, 6456–6472.
Hugonnet, R., McNabb, R., Berthier, E., Menounos, B., Nuth, C., Girod, L., Farinotti, D., Huss,
M., Dussaillant, I., Brun, F., 2021. Accelerated global glacier mass loss in the early
twenty-first century. Nature 592, 726–731.
Huss, M., 2013. Density assumptions for converting geodetic glacier volume change to mass
change. The Cryosphere 7, 877–887.
Hyndman, R., Cassidy, J., Adams, J., Rogers, G., Mazzotti, S., 2005. Earthquakes and seismic
hazard in the Yukon-Beaufort-Mackenzie. CSEG Rec. 30, 32–66.
IPCC, 2022. In: H.-O. Pörtner, D.C. Roberts, M. Tignor, E.S. Poloczanska, K. Mintenbeck, A.
Alegría, M. Craig, S. Langsdorf, S. Löschke, V. Möller, A. Okem, B. Rama (eds.),
Climate Change 2022: Impacts, Adaptation, and Vulnerability. Contribution of Working
Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate
Change. Camb. Univ. Press Press.
IPCC, 2021. In: Masson-Delmotte, V., P. Zhai, A. Pirani, S.L. Connors, C. Péan, S. Berger, N.
Caud, Y. Chen, L. Goldfarb, M.I. Gomis, M. Huang, K. Leitzell, E. Lonnoy, J.B.R.
Matthews, T.K. Maycock, T. Waterfield, O. Yelekçi, R. Yu, and B. Zhou (Eds.), Climate
Change 2021: The Physical Science Basis. Contribution of Working Group I to the Sixth
Assessment Report of the Intergovernmental Panel on Climate Change. Camb. Univ.
Press Press.
IPCC, 2013. In: Stocker, T.F., D. Qin, G.-K. Plattner, M. Tignor, S.K. Allen, J. Boschung, A.
Nauels, Y. Xia, V. Bex and P.M. Midgley (Eds.), Climate Change 2013: The Physical
Science Basis: Contribution of Working Group I to the Fifth Assessment Report of the
Intergovernmental Panel on Climate Change. Camb. Univ. Press Camb. U. K. N. Y. NY
USA.
Jackson, L., Isobe, J., 1990. Rock avalanches in the Pelly Mountains, Yukon Territory. Geol.
Surv. Can. Pap. 263–269.
Jackson, L.E., Jr., Morison, S.R., McKenna, K., 1993. Surficial geology, Prevost River, Yukon
Territory. Map 1833A.
Jacoby, G., Cook, E., 1981. Past temperature variations inferred from a 400-year tree-ring
chronology from Yukon Territory, Canada. Arct. Alp. Res. 13, 409–418.
Jóhannesson, T., Raymond, C., Waddington, E., 1989. Time–scale for adjustment of glaciers to
changes in mass balance. J. Glaciol. 35, 355–369.
Kaplan, M., Schaefer, J., Strelin, J.A., Denton, G., Anderson, R., Vandergoes, M., Finkel, R.,
Schwartz, R., Travis, S., Garcia, J., 2016. Patagonian and southern South Atlantic view of
Holocene climate. Quat. Sci. Rev. 141, 112–125.
Kaufman, D.S., Broadman, E., 2023. Revisiting the Holocene global temperature conundrum.
Nature 614, 425–435.
King, A.C., Bauska, T.K., Brook, E.J., Kalk, M., Nehrbass-Ahles, C., Wolff, E.W., Strawson, I.,
Rhodes, R.H., Osman, M.B., 2024. Reconciling ice core CO2 and land-use change
following New World-Old World contact. Nat. Commun. 15, 1735.
Kirkbride, M., Winkler, S., 2012. Correlation of Late Quaternary moraines: impact of climate
variability, glacier response, and chronological resolution. Quat. Sci. Rev. 46, 1–29.
94

Koch, A., Brierley, C., Maslin, M.M., Lewis, S.L., 2019. Earth system impacts of the European
arrival and Great Dying in the Americas after 1492. Quat. Sci. Rev. 207, 13–36.
Koch, J., Menounos, B., Clague, J.J., 2009. Glacier change in Garibaldi Provincial Park,
southern Coast Mountains, British Columbia, since the Little Ice Age. Glob. Planet.
Change 66, 161–178.
Köhler, P., Nehrbass-Ahles, C., Schmitt, J., Stocker, T.F., Fischer, H., 2017. A 156 kyr smoothed
history of the atmospheric greenhouse gases CO 2, CH 4, and N 2 O and their radiative
forcing. Earth Syst. Sci. Data 9, 363–387.
Kokelj, S.V., Lantz, T.C., Tunnicliffe, J., Segal, R., Lacelle, D., 2017. Climate-driven thaw of
permafrost preserved glacial landscapes, northwestern Canada. Geology 45, 371–374.
Lakeman, T.R., Clague, J.J., Menounos, B., 2008. Advance of alpine glaciers during final retreat
of the Cordilleran ice sheet in the Finlay River area, northern British Columbia, Canada.
Quat. Res. 69, 188–200.
Lal, D., 1991. Cosmic ray labeling of erosion surfaces: in situ nuclide production rates and
erosion models. Earth Planet. Sci. Lett. 104, 424–439.
Lal, D., Arnold, J.R., Honda, M., 1960. Cosmic-ray production rates of Be 7 in oxygen, and P
32, P 33, S 35 in argon at mountain altitudes. Phys. Rev. 118, 1626.
Lamp, J.L., Young, N.E., Koffman, T., Schimmelpfennig, I., Tuna, T., Bard, E., Schaefer, J.M.,
2019. Update on the cosmogenic in situ 14C laboratory at the Lamont-Doherty Earth
Observatory. Nucl. Instrum. Methods Phys. Res. Sect. B Beam Interact. Mater. At. 456,
157–162.
LDEO, 2012a. Separation and Purification of Quartz from Whole Rock.
LDEO, 2012b. Extraction of Beryllium from Quartz.
Leonard, E.M., 1997. The relationship between glacial activity and sediment production:
evidence from a 4450-year varve record of neoglacial sedimentation in Hector Lake,
Alberta, Canada. J. Paleolimnol. 17, 319–330.
Leys, C., Ley, C., Klein, O., Bernard, P., Licata, L., 2013. Detecting outliers: Do not use
standard deviation around the mean, use absolute deviation around the median. J. Exp.
Soc. Psychol. 49, 764–766.
Lifton, N., 2016. Implications of two Holocene time-dependent geomagnetic models for
cosmogenic nuclide production rate scaling. Earth Planet. Sci. Lett. 433, 257–268.
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.
Liu, Z., Otto-Bliesner, B., He, F., Brady, E., Tomas, R., Clark, P., Carlson, A., Lynch-Stieglitz,
J., Curry, W., Brook, E., 2009. Transient simulation of last deglaciation with a new
mechanism for Bølling-Allerød warming. science 325, 310–314.
Luckman, B., Osborn, G., 1979. Holocene glacier fluctuations in the middle Canadian Rocky
Mountains. Quat. Res. 11, 52–77.
Luckman, B.H., 2000. The little ice age in the Canadian Rockies. Geomorphology 32, 357–384.
Luckman, B.H., Holdsworth, G., Osborn, G.D., 1993. Neoglacial glacier fluctuations in the
Canadian Rockies. Quat. Res. 39, 144–153.
Lukas, S., 2011. Ice-cored moraines. Encycl. Snow Ice Glaciers 616–619.
Marshall, S.J., 2011. The cryosphere. Princeton University Press.

95

Marshall, S.J., White, E.C., Demuth, M.N., Bolch, T., Wheate, R., Menounos, B., Beedle, M.J.,
Shea, J.M., 2011. Glacier water resources on the eastern slopes of the Canadian Rocky
Mountains. Can. Water Resour. J. 36, 109–134.
Marzeion, B., Cogley, J.G., Richter, K., Parkes, D., 2014. Attribution of global glacier mass loss
to anthropogenic and natural causes. Science 345, 919–921.
Massey Jr, F.J., 1951. The Kolmogorov-Smirnov test for goodness of fit. J. Am. Stat. Assoc. 46,
68–78.
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.
McKenney, D.W., Hutchinson, M.F., Papadopol, P., Lawrence, K., Pedlar, J., Campbell, K.,
Milewska, E., Hopkinson, R.F., Price, D., Owen, T., 2011. Customized spatial climate
models for North America. Bull. Am. Meteorol. Soc. 92, 1611–1622.
Menounos, B., Clague, J.J., Clarke, G.K., Marcott, S.A., Osborn, G., Clark, P.U., Tennant, C.,
Novak, A.M., 2013a. 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., Davis, P.T., Ponce, F., Goehring, B., Maurer, M.,
Rabassa, J., Coronato, A., Marr, R., 2013b. Latest Pleistocene and Holocene glacier
fluctuations in southernmost Tierra del Fuego, Argentina. Quat. Sci. Rev. 77, 70–79.
Menounos, B., Goehring, B.M., Osborn, G., Margold, M., Ward, B., Bond, J., Clarke, G.K.,
Clague, J.J., Lakeman, T., Koch, 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., 2019. Heterogeneous changes in western North American glaciers
linked to decadal variability in zonal wind strength. Geophys. Res. Lett. 46, 200–209.
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.
Menounos, B., Pelto, B.M., Fleming, S.W., Moore, R.D., Weber, F., Hutchinson, D., Brahney, J.,
2020. Glaciers in the Canadian Columbia Basin, Technical Report.
Menounos, B., Schiefer, E., Slaymaker, O., 2006. Nested temporal suspended sediment yields,
green lake basin, British Columbia, Canada. Geomorphology 79, 114–129.
Miller, G.H., Brigham-Grette, J., Alley, R., Anderson, L., Bauch, H.A., Douglas, M., Edwards,
M., Elias, S., Finney, B., Fitzpatrick, J.J., 2010. Temperature and precipitation history of
the Arctic. Quat. Sci. Rev. 29, 1679–1715.
Miller, G.H., Geirsdóttir, Á., Zhong, Y., Larsen, D.J., Otto‐Bliesner, B.L., Holland, M.M.,
Bailey, D.A., Refsnider, K.A., Lehman, S.J., Southon, J.R., 2012. Abrupt onset of the
Little Ice Age triggered by volcanism and sustained by sea‐ice/ocean feedbacks.
Geophys. Res. Lett. 39.
Miller, G.H., Lehman, S.J., Refsnider, K.A., Southon, J.R., Zhong, Y., 2013. Unprecedented
recent summer warmth in Arctic Canada. Geophys. Res. Lett. 40, 5745–5751.
Milner, A.M., Khamis, K., Battin, T.J., Brittain, J.E., Barrand, N.E., Füreder, L., Cauvy-Fraunié,
S., Gíslason, G.M., Jacobsen, D., Hannah, D.M., 2017. Glacier shrinkage driving global
changes in downstream systems. Proc. Natl. Acad. Sci. 114, 9770–9778.
Monger, J., Price, R., 2002. The Canadian Cordillera: geology and tectonic evolution. CSEG
Rec. 27, 17–36.
96

Mood, B.J., Smith, D.J., 2015. Holocene glacier activity in the British Columbia Coast
Mountains, Canada. Quat. Sci. Rev. 128, 14–36.
Moore, R., Fleming, S., 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. Int. J. 23, 42–61.
Mukherjee, K., Menounos, B., Shea, J., Mortezapour, M., Ednie, M., Demuth, M.N., 2023.
Evaluation of surface mass-balance records using geodetic data and physically-based
modelling, Place and Peyto glaciers, western Canada. J. Glaciol. 69, 665–682.
Myhre, G., Highwood, E.J., Shine, K.P., Stordal, F., 1998. New estimates of radiative forcing
due to well mixed greenhouse gases. Geophys. Res. Lett. 25, 2715–2718.
Naegeli, K., Huss, M., Hoelzle, M., 2019. Change detection of bare-ice albedo in the Swiss Alps.
The Cryosphere 13, 397–412.
Nesje, A., Dahl, S.O., Bakke, J., 2004. Were abrupt Lateglacial and early-Holocene climatic
changes in northwest Europe linked to freshwater outbursts to the North Atlantic and
Arctic Oceans? The Holocene 14, 299–310.
Nesje, A., Matthews, J.A., Dahl, S.O., Berrisford, M.S., Andersson, C., 2001. Holocene glacier
fluctuations of Flatebreen and winter-precipitation changes in the Jostedalsbreen region,
western Norway, based on glaciolacustrine sediment records. The Holocene 11, 267–280.
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. Sect. B
Beam Interact. Mater. At. 258, 403–413.
Nuth, C., Kääb, A., 2011. Co-registration and bias corrections of satellite elevation data sets for
quantifying glacier thickness change. The Cryosphere 5, 271–290.
Nye, J.F., 1960. The response of glaciers and ice-sheets to seasonal and climatic changes. Proc.
R. Soc. Lond. Ser. Math. Phys. Sci. 256, 559–584.
Obu, J., Westermann, S., Bartsch, A., Berdnikov, N., Christiansen, H.H., Dashtseren, A.,
Delaloye, R., Elberling, B., Etzelmüller, B., Kholodov, A., 2019. Northern Hemisphere
permafrost map based on TTOP modelling for 2000–2016 at 1 km2 scale. Earth-Sci. Rev.
193, 299–316.
Oerlemans, J., 1994. Quantifying global warming from the retreat of glaciers. Science 264, 243–
245.
Osborn, G., McCarthy, D., LaBrie, A., Burke, R., 2015. Lichenometric dating: science or
pseudo-science? Quat. Res. 83, 1–12.
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., Robinson, B., Luckman, B., 2001. Holocene and latest Pleistocene fluctuations of
Stutfield Glacier, Canadian Rockies. Can. J. Earth Sci. 38, 1141–1155.
Osterberg, E.C., Mayewski, P.A., Fisher, D.A., Kreutz, K.J., Maasch, K.A., Sneed, S.B., Kelsey,
E., 2014. Mount Logan ice core record of tropical and solar influences on Aleutian Low
variability: 500–1998 AD. J. Geophys. Res. Atmospheres 119, 11–189.
Owens, M.J., Lockwood, M., Hawkins, E., Usoskin, I., Jones, G.S., Barnard, L., Schurer, A.,
Fasullo, J., 2017. The Maunder minimum and the Little Ice Age: an update from recent
reconstructions and climate simulations. J. Space Weather Space Clim. 7, A33.

97

Pendleton, S.L., Briner, J.P., Kaufman, D.S., Zimmerman, S.R., 2017. Using cosmogenic 10Be
exposure dating and lichenometry to constrain Holocene glaciation in the central Brooks
Range, Alaska. Arct. Antarct. Alp. Res. 49, 115–132.
Phillips, F.M., Argento, D.C., Balco, G., Caffee, M.W., Clem, J., Dunai, T.J., Finkel, R.,
Goehring, B., Gosse, J.C., Hudson, A.M., 2016. The CRONUS-Earth project: a synthesis.
Quat. Geochronol. 31, 119–154.
Porter, C., Howat, I., Noh, M.-J., Husby, E., Khuvis, S., Danish, E., Tomko, K., Gardiner, J.,
Negrete, A., Yadav, B., Klassen, J., Kelleher, C., Cloutier, M., Bakker, J., Enos, J.,
Arnold, G., Bauer, G., Morin, P., 2023. ArcticDEM - Mosaics, Version 4.1.
https://doi.org/10.7910/DVN/3VDC4W
Porter, C., Howat, I., Noh, M.-J., Husby, E., Khuvis, S., Danish, E., Tomko, K., Gardiner, J.,
Negrete, A., Yadav, B., Klassen, J., Kelleher, C., Cloutier, M., Bakker, J., Enos, J.,
Arnold, G., Bauer, G., Morin, P., 2022. ArcticDEM - Strips, Version 4.1.
https://doi.org/10.7910/DVN/C98DVS
Porter, S.C., 1986. Pattern and forcing of Northern Hemisphere glacier variations during the last
millennium. Quat. Res. 26, 27–48.
Porter, S.C., Denton, G.H., 1967. Chronology of neoglaciation in the North American Cordillera.
Am. J. Sci. 265, 177–210.
Rasmussen, S.O., Andersen, K.K., Svensson, A., Steffensen, J.P., Vinther, B.M., Clausen, H.B.,
Siggaard‐Andersen, M., Johnsen, S.J., Larsen, L.B., Dahl‐Jensen, D., 2006. A new
Greenland ice core chronology for the last glacial termination. J. Geophys. Res.
Atmospheres 111.
Reichert, B.K., Bengtsson, L., Oerlemans, J., 2002. Recent glacier retreat exceeds internal
variability. J. Clim. 15, 3069–3081.
Renssen, H., Seppä, H., Heiri, O., Roche, D.M., Goosse, H., Fichefet, T., 2009. The spatial and
temporal complexity of the Holocene thermal maximum. Nat. Geosci. 2, 411–414.
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., 2006. Expansion of alpine glaciers in Pacific North America
in the first millennium AD. Geology 34, 57–60.
Reza, A.M., 2004. Realization of the contrast limited adaptive histogram equalization (CLAHE)
for real-time image enhancement. J. VLSI Signal Process. Syst. Signal Image Video
Technol. 38, 35–44.
Riahi, K., Van Vuuren, D.P., Kriegler, E., Edmonds, J., O’neill, B.C., Fujimori, S., Bauer, N.,
Calvin, K., Dellink, R., Fricko, O., 2017. The Shared Socioeconomic Pathways and their
energy, land use, and greenhouse gas emissions implications: An overview. Glob.
Environ. Change 42, 153–168.
Rinterknecht, V.R., 2003. Cosmogenic 10 Be chronology for the last deglaciation of the southern
Scandinavian Ice Sheet.
Roe, G.H., O’Neal, M.A., 2009. The response of glaciers to intrinsic climate variability:
observations and models of late-Holocene variations in the Pacific Northwest. J. Glaciol.
55, 839–854.
Rounce, D.R., Hock, R., Maussion, F., Hugonnet, R., Kochtitzky, W., Huss, M., Berthier, E.,
Brinkerhoff, D., Compagno, L., Copland, L., 2023. Global glacier change in the 21st
century: Every increase in temperature matters. Science 379, 78–83.
Rousseeuw, P.J., Croux, C., 1993. Alternatives to the median absolute deviation. J. Am. Stat.
Assoc. 88, 1273–1283.
98

Ruddiman, W.F., Fuller, D.Q., Kutzbach, J.E., Tzedakis, P.C., Kaplan, J.O., Ellis, E.C., Vavrus,
S.J., Roberts, C.N., Fyfe, R., He, F., 2016. Late Holocene climate: Natural or
anthropogenic? Rev. Geophys. 54, 93–118.
Ryder, J., Thomson, B., 1986. Neoglaciation in the southern Coast Mountains of British
Columbia: chronology prior to the late Neoglacial maximum. Can. J. Earth Sci. 23, 273–
287.
Schaefer, J.M., Codilean, A.T., Willenbring, J.K., Lu, Z.-T., Keisling, B., Fülöp, R.-H., Val, P.,
2022. Cosmogenic nuclide techniques. Nat. Rev. Methods Primer 2, 1–22.
Schaefer, J.M., Denton, G.H., Kaplan, M., Putnam, A., Finkel, R.C., Barrell, D.J., Andersen,
B.G., Schwartz, R., Mackintosh, A., Chinn, T., 2009. High-frequency Holocene glacier
fluctuations in New Zealand differ from the northern signature. science 324, 622–625.
Schimmelpfennig, I., Schaefer, J.M., Lamp, J., Godard, V., Schwartz, R., Bard, E., Tuna, T.,
Akçar, N., Schlüchter, C., Zimmerman, S., 2022. Glacier response to Holocene warmth
inferred from in situ 10 Be and 14 C bedrock analyses in Steingletscher’s forefield
(central Swiss Alps). Clim. Past 18, 23–44.
Segal, R.A., Lantz, T.C., Kokelj, S.V., 2016. Acceleration of thaw slump activity in glaciated
landscapes of the Western Canadian Arctic. Environ. Res. Lett. 11, 034025.
Sigl, M., Toohey, M., McConnell, J.R., Cole-Dai, J., Severi, M., 2022. Volcanic stratospheric
sulfur injections and aerosol optical depth during the Holocene (past 11 500 years) from a
bipolar ice-core array. Earth Syst. Sci. Data 14, 3167–3196.
Sigl, M., Toohey, M., McConnell, J.R., Cole-Dai, J., Severi, M., 2021. HolVol: Reconstructed
volcanic stratospheric sulfur injections and aerosol optical depth for the Holocene (9500
BCE to 1900 CE).
Smith, C.A., Bennett, B.A. (Eds.), 2004. Ecoregions of the Yukon Territory: Biophysical
properties of Yukon landscapes. Summerland, BC: Agriculture and Agri-food Canada.
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., 2015. Holocene glacier fluctuations.
Quat. Sci. Rev. 111, 9–34.
Steinhilber, F., Beer, J., Fröhlich, C., 2009. Total solar irradiance during the Holocene. Geophys.
Res. Lett. 36.
Stroeven, A.P., Fabel, D., Codilean, A.T., Kleman, J., Clague, J.J., Miguens-Rodriguez, M., Xu,
S., 2010. Investigating the glacial history of the northern sector of the Cordilleran Ice
Sheet with cosmogenic 10Be concentrations in quartz. Quat. Sci. Rev. 29, 3630–3643.
Stroeven, A.P., Fabel, D., Margold, M., Clague, J.J., Xu, S., 2014. Investigating absolute
chronologies of glacial advances in the NW sector of the Cordilleran Ice Sheet with
terrestrial in situ cosmogenic nuclides. Quat. Sci. Rev. 92, 429–443.
Tam, B.Y., Szeto, K., Bonsal, B., Flato, G., Cannon, A.J., Rong, R., 2019. CMIP5 drought
projections in Canada based on the Standardized Precipitation Evapotranspiration Index.
Can. Water Resour. JournalRevue Can. Ressour. Hydr. 44, 90–107.
Tennant, C., Menounos, B., 2013. Glacier change of the Columbia Icefield, Canadian Rocky
Mountains, 1919–2009. J. Glaciol. 59, 671–686.
Thiagarajan, N., Subhas, A.V., Southon, J.R., Eiler, J.M., Adkins, J.F., 2014. Abrupt preBølling–Allerød warming and circulation changes in the deep ocean. Nature 511, 75–78.
Tierney, J.E., Poulsen, C.J., Montañez, I.P., Bhattacharya, T., Feng, R., Ford, H.L., Hönisch, B.,
Inglis, G.N., Petersen, S.V., Sagoo, N., 2020. Past climates inform our future. science
370, eaay3701.
99

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.
Tomkins, M.D., Dortch, J.M., Hughes, P.D., Huck, J.J., Pallàs, R., Rodés, Á., Allard, J.L.,
Stimson, A.G., Bourlès, D., Rinterknecht, V., 2021. Moraine crest or slope: An analysis
of the effects of boulder position on cosmogenic exposure age. Earth Planet. Sci. Lett.
570, 117092.
Tomkins, M.D., Huck, J.J., Dortch, J.M., Hughes, P.D., Kirbride, M.P., Barr, I.D., 2018. Schmidt
Hammer exposure dating (SHED): Calibration procedures, new exposure age data and an
online calculator. Quat. Geochronol. 44, 55–62.
Toohey, M., Sigl, M., 2019. Reconstructed volcanic stratospheric sulfur injections and aerosol
optical depth, 500 BCE to 1900 CE, version 3.
https://doi.org/10.26050/WDCC/eVolv2k_v3
Urbanski, J.A., 2018. A GIS tool for two-dimensional glacier-terminus change tracking. Comput.
Geosci. 111, 97–104.
van Dijk, E.J., Jungclaus, J., Sigl, M., Timmreck, C., Krüger, K., 2024. High-frequency climate
forcing causes prolonged cold periods in the Holocene. Commun. Earth Environ. 5, 242.
Walker, A., Hidy, A.J., Zimmerman, S., Tumey, S.J., Brown, T.A., Braumann, S.M., Schwartz,
R., Schaefer, J.M., 2022. Improved precision for 10Be-based chronology by optimizing
carrier addition.
Wanner, H., Beer, J., Bütikofer, J., Crowley, T.J., Cubasch, U., Flückiger, J., Goosse, H.,
Grosjean, M., Joos, F., Kaplan, J.O., 2008. Mid-to Late Holocene climate change: an
overview. Quat. Sci. Rev. 27, 1791–1828.
Wanner, H., Solomina, O., Grosjean, M., Ritz, S.P., Jetel, M., 2011. Structure and origin of
Holocene cold events. Quat. Sci. Rev. 30, 3109–3123.
Wiles, G.C., D’Arrigo, R.D., Villalba, R., Calkin, P.E., Barclay, D.J., 2004. Century‐scale solar
variability and Alaskan temperature change over the past millennium. Geophys. Res.
Lett. 31.
Young, N.E., Briner, J.P., Kaufman, D.S., 2009. Late Pleistocene and Holocene glaciation of the
Fish Lake valley, northeastern Alaska Range, Alaska. J. Quat. Sci. Publ. Quat. Res.
Assoc. 24, 677–689.
Zemp, M., Hoelzle, M., Haeberli, W., 2009. Six decades of glacier mass-balance observations: a
review of the worldwide monitoring network. Ann. Glaciol. 50, 101–111.

100