# Olivine Time-Capsules Constrain the Pre-Eruptive History of Holocene Basalts, Mount Meager Volcanic Complex, British Columbia, Canada

S.M. AUFRÈRE<fnref n="1" />, G. WILLIAMS-JONES<fnref n="1" />, S. MOUNE<fnref n="2" />, D.J. MORGAN<fnref n="3" />, N. VIGOUROUX<fnref n="1" />,<fnref n="4" /> and J.K. RUSSELL<fnref n="5" />

The Canadian segment of the Cascade Volcanic Arc (i.e. the Garibaldi Volcanic Belt) comprises more than 100 eruptive centres, spanning the entire Quaternary period (Pleistocene to Holocene in age), and with deposits ranging in composition from alkaline basalt to rhyolite. At least one of the volcanoes is currently active; Mount Meager / *Quelwěwelüsten* erupted explosively 2360 years BP and has ongoing fumarolic activity. Long-term forecasting of eruption frequency and style depends on reconstruction of the history and timescales of magmatic processes preceding previous volcanic eruptions. Utilising diffusion chronometry, we investigate the Mount Meager Volcanic Complex focusing on Holocene olivine-phyric basalts (Lillooet Glacier basalts) exposed by the retreat of the Lillooet Glacier. We identify two distinct olivine populations in samples of quenched, glassy basalt lavas that record different magmatic processes and histories. Glomerocrysts of Fo₈₃ olivine phenocrysts, entrained and transported by a hot mafic input, form Population 1. These exhibit resorption and normally zoned outermost rim compositions of Fo₇₆-7₈; a third of them also show interior reverse compositional zoning. A second population of skeletal microphenocrysts have the same composition as the phenocryst rims (i.e. Fo₇₆-7₈) and are in equilibrium with the adjacent matrix glass. We estimate the pre-eruptive temperature-fO₂ conditions in a shallow reservoir (100 MPa; ~3 km) for a melt with H₂O content of 0.5 to 1 wt % as ~1097°C to 1106°C (± 30°C), and NNO + 0.5 (±1.1), respectively. Using these input parameters, we report Fe-Mg diffusion chronometry results for 234 normally zoned profiles from 81 olivine phenocrysts. Diffusion modelling of compositional profiles in oriented crystals indicates pre-eruptive magmatic residence times of 1 to 3 months. These remarkably short residence times in shallow reservoirs prior to eruption suggest very short periods of unrest may precede future eruptions.

**Key words:** Basalt; Cascade Volcanic Arc; Diffusion chronometry; Mount Meager Volcanic Complex; Olivine

## INTRODUCTION

Reconstructing the history of magmatic processes that preceded volcanic unrest and eruption is crucial in forecasting and understanding future eruptions. Many arc magmas show evidence for multiple pre-eruptive events and processes that could include magma storage and crystallisation at shallow crustal levels, influx of more primitive magma, magma mixing, or the assimilation of host crustal rocks. However, it is often challenging to identify which magmatic process(es) primed the system and, ultimately, triggered an eruption (e.g. Ubide et al., 2019; Mourey et al., 2023).

Here we study olivine crystals representative of the Lillooet Glacier (LG) basalts in the Garibaldi Volcanic Belt (GVB), i.e. the northern segment of the Cascade Volcanic Arc (Fig. 1a, b). These basalts, located within the Mount Meager Volcanic Complex (MMVC), are among the youngest deposits in the area and likely Holocene in age; Wilson & Russell (2017) found they were too young to date by ⁴⁰Ar/³⁹Ar geochronology. Our investigation commences with a concise overview of the various mineral phases

observed in the LG basalts, followed by a detailed analysis of the two primary olivine populations observed: resorbed phenocrysts with Fo₈₃ core compositions, and skeletal microphenocrysts with Fo₇₆-7₈ core compositions. The compositional zoning patterns in olivine phenocrysts are also quantified and used for Fe-Mg diffusion chronometry. Modelling of these data informs on the magma storage conditions and residence times immediately preceding eruption and are relevant to the timescales of volcanic unrest to be expected prior to future eruptions.

This approach of estimating timescales between chemical disequilibrium and eruption has been shown to agree with shallow and deep seismicity, surface deformation, thermal anomalies and gas ratio changes (e.g. Kahl et al., 2011, 2013, 2022, 2023a; Saunders et al., 2012; Pankhurst et al., 2018; Rasmussen et al., 2018; Mourey et al., 2023; Voloschina et al., 2023). Thus, diffusion chronometry presents itself as a promising tool for estimating perturbation-to-eruption timescales in volcanoes lacking monitoring systems, potentially aiding in the understanding and anticipation of future volcanic activity. To date, diffusion chronometry has been applied to a few volcanic centres in the Cascade Volcanic

<footnote n="1">¹Centre for Natural Hazards Research, Department of Earth Sciences, TASC 1, Simon Fraser University, 8888 University Drive, Burnaby, BC V5A 1S6, Canada</footnote>

<footnote n="2">²Laboratoire Magmas et Volcans, CNRS-OPGC-IRD, Université Clermont Auvergne, Campus Universitaire des Cézeaux, 6 Avenue Blaise Pascal, Aubière 63170, France</footnote>

<footnote n="3">³School of Earth and Environment, University of Leeds, Woodhouse, Leeds LS2 9JT, United Kingdom</footnote>

<footnote n="4">⁴Department of Earth and Environmental Sciences, Douglas College, 700 Royal Avenue, New Westminster, BC V3M 5Z5, Canada</footnote>

<footnote n="5">⁵Department of Earth, Ocean and Atmospheric Sciences, Volcanology and Petrology Laboratory, University of British Columbia, 2207 Main Mall #2020, Vancouver, BC V6T 1Z4, Canada</footnote>

<footnote>*Corresponding author. E-mail: sarah_aufrere@sfu.ca.</footnote>
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Geological setting of the Mount Meager Volcanic Complex (MMVC: larger triangle in inset). a) Map of the Cascades Volcanic Arc (CVA), that extends from northern California, USA, to southwestern British Columbia, Canada (modified from Harris et al., 2023; Harris & Russell, 2022; Wilson & Russell, 2017). The northern segment of the CVA is formed by the Garibaldi Volcanic Belt (GVB). b) Geological map of the MMVC based on Harris et al. (2022, 2023) showing the extent of the MMVC volcanic deposits. These are split into the mafic, intermediate, felsic or undivided felsic-to-intermediate compositions. When several ages were available for one location or assemblage, the youngest (given this age has a small standard deviation) is displayed in italic below the name of the assemblage (Devastator Creek, Eastern and Western Lillooet Ridge, Elaho Valley, Lillooet River, Mosaic Ridge: Harris et al., 2023; Cracked Mountain: Harris et al., 2023; Read, 1990; Lillooet Glacier: Wilson & Russell, 2017, 2018; Pebble Creek Formation: Hickson et al., 1999; Read, 1990; Plinth Block and Ash: Russell et al., 2021). Location of the Lillooet Glacier AW-15-070 sample shown by a whitestar, and location of the Job Glacier fumaroles as white flowing lines. Abbreviations: Expl, Explorer Plate; BR, Bridge River; SG, Salal Glacier; MC, Mount Cayley; MG, Mount Garibaldi; MB, Mount Baker; GP, Glacier Peak; MR, Mount Rainier; MSH, Mount St. Helens; MA, Mount Adams; MH, Mount Hood; MJ, Mount Jefferson; TS, Three Sisters; N, Newberry; CL, Crater Lake; MS, Mount Shasta; LP, Lassen Peak.</figcaption></figure>

Arc: Mount St. Helens (Li in plagioclase, Kent et al., 2007; Fe-Mg in orthopyroxene, Saunders et al., 2012), Mount Hood (Mg in plagioclase; Kent et al., 2010), Blue Lake Crater (Fe-Mg and Ni in olivine; Johnson & Cashman, 2020), Mount Shasta (Fe-Mg in clinopyroxene and orthopyroxene; Phillips & Till, 2022), Box Canyon (Ni in olivine; Hollyday et al., 2020) and Cinder Cone (Fe-Mg and Ni in olivine; Walowski et al., 2019). These studies revealed mixing-to-eruption timescales ranging from weeks to decades (Kent et al., 2010;

Saunders et al., 2012; Hollyday et al., 2020; Johnson & Cashman, 2020; Phillips & Till, 2022), and residence time in an evolving magma from weeks to years (Walowski et al., 2019).

In this study, the term *phenocryst* refers to its original non-genetic definition, i.e. crystals large enough to be seen with the unaided eye (≥ 30 µm in width) (Iddings, 1892). Likewise, the term *microphenocryst* is used to denote phenocrysts <500 µm in length, following Zellmer's (2021) classification. The Zellmer (2021 and
references therein) genetic terminology is used when discussing the origin of the different crystal populations observed here. An *autocryst* refers to a crystal (or part of one) that nucleated from the last carrier melt and is, thus, in equilibrium with the groundmass or glass. An *antecryst* describes an older crystal that originated from the same magmatic system and was subsequently incorporated into the carrier melt; it is thus chemically similar but likely shows disequilibrium textures. Another term, *philocryst* (Mangler et al., 2020), is used to refer to all antecrysts remobilised from various crystal mushes throughout a transcrustal magmatic system (Cashman et al., 2017), typical of subduction zones. Finally, *xenocryst* denotes a crystal foreign to the magmatic system (i.e. sourced from a wall or basement rock of markedly different composition) that was incorporated and now may show a reaction core and/or an overgrowth rim, the latter being stable in the magma where it currently resides (Jerram & Martin, 2008; Zellmer, 2021).

# GEOLOGICAL CONTEXT

## The Garibaldi Volcanic Belt

The southernmost volcanoes of British Columbia result from the subduction of the Juan de Fuca plate beneath the North American plate (Fig. 1a; du Bray & John, 2011; Mullen et al., 2017). They form the GVB, the northern segment of the Cascade Volcanic Arc, which extends down to NW Washington state (Fig. 1a). The Canadian portion of the GVB consists of >100 eruptive centres including three major long-lived stratovolcano complexes, Nch'kay / Mount Garibaldi, Tak'takmu'yin t'l'a in7in'axa7en / Mount Cayley, and *Qwel'qwelústen* / Mount Meager, with compositions ranging from alkaline basalts to rhyolite (e.g. Russell et al., 2023). There are notably lower eruption rates and erupted volumes compared to the southern Cascades, which most likely reflects a discrepancy in the underlying subduction system (Mullen et al., 2017). Indeed, trace element and isotope ratios show greater mantle depletion and higher slab-derived sediment contributions to the south (Mullen & Weis, 2015). Such north-south variations have been attributed to the Nootka fault zone, at the northern end of the GVB (Mullen & Weis, 2015; Venugopal et al., 2020), along which the Explorer plate detached from the Juan de Fuca plate at ~4 Ma (Riddihough, 1984). There, flows of enriched (OIB-like) asthenospheric mantle, coming from beneath the slab and triggered by slab rollback, are drawn into the mantle wedge, and become diluted southwards (Thorkelson et al., 2011; Mullen & Weis, 2015). However, recent geophysical studies suggest that the Nootka fault zone may now follow a more northerly trajectory than initially thought (Merrill et al., 2022).

Regional mapping of the GVB volcanoes and dating is still quite recent (e.g. Kelman, 2005; Wilson & Russell, 2018; Harris et al., 2020; Russell et al., 2021), although valuable work dates back to 1958 (Mathews, 1958a, 1958b) and the late 1970s (e.g. Easterbrook, 1975; Green, 1977; Read, 1977b). Nevertheless, only a handful of petrological studies have been conducted throughout the different GVB volcanic fields (Mullen & Weis, 2013, 2015; Mullen et al., 2017; Venugopal et al., 2020) and at specific volcanic centres (Wilson & Russell, 2017; Harris & Russell, 2022). Hence, it is imperative to continue to refine our understanding of magma storage conditions (including pressure, H₂O content, and temperature) and magmatic processes (such as potential mixing events and timescales of transport, storage, and eruption) within the GVB.

## Mount Meager Volcanic Complex

The MMVC is an active, glacier-clad calc-alkaline volcanic complex located 150 km north of Vancouver in British Columbia, Canada, and forms part of the Garibaldi Volcanic Belt (Fig. 1a;

Russell et al., 2023). We define the MMVC to comprise the Mount Meager massif (*Qwel'qwelústen* in the local Lílwat First Nation language) and the surrounding volcanic centres located within a 20-km radius (Fig. 1b).

Due to long-lived activity (>1.9 Ma) and glacier retreat, Mount Meager is a highly unstable massif, with numerous large landslides originating from its flanks, including the largest in Canadian history which occurred in 2010 (e.g. Roberti et al., 2018; Connelly et al., 2024). Moreover, anecdotal reports of sulphurous smells suggest the presence of fumaroles in the valley near Job Glacier (Fig. 1b), going back decades, possibly as far back as the 1930s (M. Kelman, pers. comm., 2024). Since 2015, such fumarolic degassing can be observed at Job Glacier through glaciovolcanic caves formed largely due to climate warming and associated glacier thinning (Roberti et al., 2018; Unnsteinsson et al., 2024). This fumarolic activity contributes to Mount Meager being classified as one of two Very High threat volcanoes in Canada, comparable to Mount Baker, south of the Canada-USA border (Fig. 1a; Kelman & Wilson, 2024).

The MMVC comprises glaciovolcanic and non-glaciovolcanic deposits (Wilson & Russell, 2018), with the massif itself composed of intermediate to felsic units (green in Fig. 1b), taking their names from the nearest peak: the *Devastator Assemblage* (1900 ka), the *Pylon Assemblage* (900 ka), the *Capricorn Assemblage* (100 ka), the *Job Assemblage*, and *Plinth Peak* (57 to 200 ka) (Read, 1977b, 1990; Green et al., 1988; Wilson & Russell, 2018). The most recent eruption dated at 2360 ¹⁴C years BP, consisted of a VEI 4 sub-Plinian eruption which produced pyroclastic density currents followed by a Vulcanian episode of dacitic lava and block and ash flows (Stasiuk et al., 1996; Hickson et al., 1999; Michol et al., 2008; Andrews et al., 2014; Russell et al., 2021) leading to the emplacement of the *Pebble Creek Formation* (Fig. 1b; Clague et al., 1995; Hickson et al., 1999; Read, 1990). In addition, the Mount Meager massif is surrounded by small (< 1 km³) volcanic centres (Fig. 1b) to the northwest (*Lillooet Glacier*, Holocene; Wilson & Russell, 2017, 2018), north (*Western and Eastern Lillooet Ridge*, 1120 and 201 ka, respectively; and *Lillooet River*, 441 ka; Harris et al., 2023), east (*Plinth Block and Ash*, 24.3 ka; Russell et al., 2021), southwest (*Elaho Valley*, 106 ka; Harris et al., 2023; *Cracked Mountain*, 401 ka; Harris et al., 2023; Read, 1990; and *Devastator Creek*, 499 ka; Harris et al., 2023) and west (*Mosaic Ridge*, 115 ka; Harris et al., 2023).

While valuable work has been conducted on the MMVC since the late 1970s (Read, 1977a, 1977b), little is known about its magmatic plumbing system (Russell et al., 2023). Some recent geophysical studies (based on magnetotelluric imaging) revealed that the deeper magmatic system (~5–15 km below sea level) may have a minimum volume of 2 × 10¹² m³, with a melt fraction of ~18–32% of dacitic-to-trachytic composition at ~800°C to 900°C (Hanneson & Unsworth, 2023). Petrological studies to date include one olivine-hosted melt inclusion study conducted on *Mosaic Ridge* (Venugopal et al., 2020), which revealed polybaric olivine crystallisation (100–400 MPa) occurring at temperatures between 1057°C and 1142°C (± 50°C) and under NNO + 0.74 conditions. This investigation also identified exceptionally high H₂O levels (1.3–3.1 wt %) and CO₂ concentrations (>1500 ppm), marking them as the highest recorded across the GVB eruptions. Other studies have primarily explored the glaciovolcanic origins of specific mafic centres within the MMVC (*Cracked Mountain*, *Lillooet Glacier*; Fig. 1b; Harris & Russell, 2022; Wilson & Russell, 2017).

## Lillooet Glacier

The LG unit, located 20 km northwest of the Mount Meager massif (Fig. 1b), was recently exposed due to the ~5 km-up-valley
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Total alkali silica diagram (Le Bas *et al.*, 1986) of the whole-rock compositions from the MMVC mafic centres: Lillooet Glacier, Lillooet Ridge, Lillooet River, Cracked Mountain, Elaho Valley, Mosaic Ridge (see Fig. 1b). Whole-rock analyses are from Harris *et al.* (2023) and Wilson & Russell (2017). The white star represents the Lillooet Glacier sample (AW-15-070) used in this study.</figcaption></figure>

retreat of the Lillooet Glacier (Walker, 2003; Reyes & Clague, 2004) since the end of the Little Ice Age (~1900 CE; Menounos *et al.*, 2009). Wilson & Russell (2017) mapped and sampled the LG unit (sample AW-15-070 used here; Fig. 1b, 2). The LG basalts were shown to be glaciovolcanic in origin and likely to relate to the waning of the last (Fraser) glaciation (i.e. < 11 ka; e.g. Clague & Ward, 2011), making it the youngest mafic centre in the MMVC. Wilson & Russell's (2017) petrochemical analyses show the LG lavas to be comagmatic; thus, they share a similar composition and mineralogy, with all apparent chemical variations explained by crystal-sorting of the observed olivine + plagioclase phenocryst assemblage. A thorough description of the lithofacies, petrography and geochemistry of the LG volcanics is presented in Wilson & Russell (2017). They also constrained the LG pre-eruptive storage conditions at < 200 MPa, with an H₂O content of 0.5 to 1 wt % (Wilson & Russell, 2017). These were inferred both from thermodynamic modelling (using rhyolite-MELTS; Gualda *et al.*, 2012; Ghiorso & Gualda, 2015) involving equilibrium crystallisation of plagioclase + olivine + augite gabbroic material during storage (as seen through a cognate, gabbroic inclusion in their sample 59; see Fig. 7C-D in Wilson & Russell, 2017) and isobaric (pressure fixed at 100 MPa) fractional crystallisation of olivine and plagioclase, where H₂O contents were dictated by plagioclase saturation only (Wilson & Russell, 2017). Additionally, further simulations based on their most evolved measured glass composition lead Wilson & Russell (2017) to estimate the eruption temperature to be ~1045°C. As such, Wilson & Russell's (2017) comprehensive mapping, petrography and geochemistry of the LG lavas and their pre-eruptive magmatic condition modelling serve as an excellent foundation to the current study, which builds upon it by further constraining the LG pre-eruptive magmatic system and its timescales through an olivine lens.

<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Chemical bimodality of the two main olivine populations (Groups 1 and 2) observed within the Lillooet Glacier lavas. These are plotted against their Fo and NiO content, and Kernel density estimations plots are represented on the opposite axes. Data from Wilson & Russell (2017) and this study.</figcaption></figure>

# SAMPLING AND ANALYTICAL METHODS

## Sample selection rationale

The LG eruptive centre is remote, and access requires helicopter support. We, therefore, accessed the sample archive of Wilson & Russell (2017) for our more detailed study of the LG lavas. Samples gathered in Wilson & Russell (2017) include pillow lavas (both the rim and the core of the pillows), glassy pillow rinds, pillow breccia, volcaniclastic sand (inter-pillow material), lava tubes, and dense basalt. For diffusion chronometry, samples need to be quenched upon eruption to ensure little to no post-eruption crystal growth and diffusion; thus, using glassy pillow rinds was the most appropriate option. Hence, sample AW-15-070, the glassiest sample in the unit (see Figs. A1 and A3 in Appendix A) and consisting of glassy pillow rinds from multiple pillows, was chosen (sample collected on June 26, 2015 by A. Wilson; GPS location: 50.7471, -123.7255; elevation: 1150 m asl). In addition, glass analyses, key to geothermometry, are thus possible. The olivine grains present in sample AW-15-070 show the same textures (Fig. A3) and chemistry (see Fig. 3) as all LG samples studied by Wilson & Russell (2017) and are thus considered representative of their respective populations throughout the LG unit. Wilson & Russell's (2017) AW-15-070 thin section was used in this study and part of this sample was crushed, with 83 olivine grains handpicked, mounted together, and then polished (see Appendix B for details).

# Electron Backscatter Diffraction (EBSD)

Fe-Mg interdiffusion in olivine is anisotropic, being four to eight times faster along the c-axis than along the a- and b-axes in the temperature range 800°C to 1200°C (Dohmen *et al.*, 2007; Dohmen & Chakraborty, 2007a, 2007b). Therefore, it is essential to determine grain orientations to calculate accurate timescales. A low-resolution montage map of the grain-mount and grain orientations was obtained using the CamScan X500FE Crystal Probe Scanning Electron Microscope (SEM) at Géosciences Montpellier (CNRS, Université de Montpellier, France) equipped with an Oxford/Nordlys Electron Backscatter Diffraction (EBSD) detector.
The grain-mount montage map was used as the reference for navigation and grain identification in the following analyses (Fig. C1 in Appendix).

For grain orientation purposes, the grain-mount was not carbon-coated to ensure that the diffraction (Kikuchi) pattern from the olivine crystal lattices was not obscured. The SEM was run in low vacuum, with an accelerating voltage of 15 to 20 kV and a working distance of 20 mm. A 70° tilt, necessary to obtain an optimised diffracted signal and, therefore, better and faster analyses (Oxford Instruments, 2023), was ensured by the column being pre-tilted rather than using a tilted sample holder.

EBSD patterns were indexed automatically using the Channel5 software package from Oxford Instruments HKL. Importantly, the three Euler angles were obtained for each grain. Most EBSD software (like Channel5) use the Euler angles to define the orientation of the mineral as a sequence of rotations about specific axes. Once these Euler angles were obtained, the rotations were transformed into the resultant *a*, *b* and *c* crystal axis orientations using an in-house Excel spreadsheet (Euler-proc Excel spreadsheet; D. J. Morgan, pers. comm., 2022). The MATLAB-supported MTEX software (Bachmann et al., 2010) was then used to ensure that the orientations were correct. The MTEX code considers as an input the EBSD data, i.e. the three Euler angles following the Bunge convention, and the *x* and *y* coordinates of each grain. To keep the reference frame consistent, the *x*-axis was plotted eastwards and the *y*-axis southwards. Oriented olivine crystals were then plotted on top of the grain-mount crystals and found to match those given by the Euler-proc spreadsheet (Appendix Fig. C2).

## Electron Probe Micro-Analysis (EPMA)

Figure D1 (see Appendix D) summarises all the procedures conducted with an Electron Probe Microanalyzer (EPMA). This study utilised two EPMA instruments: a JEOL JXA-iHP200F EPMA located at the University of British Columbia (UBC) in Vancouver, Canada, and a CAMECA SXFive-TACTIS EPMA situated at the Laboratoire Magmas et Volcans (LMV) in Clermont-Ferrand, France.

Point analyses of the glass and all the mineral phases observed in the LG basalts (Fig. A1) i.e. olivine, spinel, clinopyroxene and plagioclase were carried out on the thin section and olivine grain-mount. All the analytical conditions and errors can be found in Appendix D, and the resulting analyses in this study's complementary dataset (Aufrère & Williams-Jones, 2024). In the grain-mount, the core of each olivine grain was analysed, as well as some associated spinel inclusions. When preserved during the crushing process, the glass in direct contact with the olivine rims (analysed as part of the profile end-member, see Compositional profiles section) was also analysed, so that the pair could be used for geothermometry. Likewise, the core and rim of several olivine microphenocrysts, as well as surrounding glass were analysed throughout the thin section. Glass analyses were also carried out on the large, apparent melt inclusions (i.e. features that resemble melt inclusions in the studied sections but may, in fact, be embayments) and on embayments (where connectivity to the surrounding melt is visible), to determine whether they had the same composition as the groundmass glass.

A key component of this study is to determine whether the embayed texture of the olivine grains was related to rapid growth or dissolution. Although the olivine growth record is hidden by Fe-Mg zoning in BSE images, it can be revealed by phosphorus (P) zoning (Welsch et al., 2014). P-zoning in rapidly grown, skeletal olivine crystals would reveal feathery, primary dendrite branches and the very beginning of the infilling stage (the olivine texture can be referred to as 'spongy'; e.g. Fig.1.2 in Welsch et al., 2014). In

contrast, resorbed olivine would show those branches irregularly crosscut (e.g. Fig. 8 in Milman-Barris et al., 2008) and/or no correlation between the crystal edges and the zones of P-enrichment, and, in the case of several dissolution events, internal dissolution fronts (e.g. Fig. 2 in Mourey et al., 2022). For these reasons, X-ray maps of Fe and P of selected olivine, displaying embayments or cuspate margins, were collected using the UBC EPMA with an accelerating voltage of 20 kV, a probe current of 800 nA, and dwell time per pixel of 240 ms. The beam was fully focused, which at these conditions resulted in a physical size of 250 nm (A. Von Der Handt, pers. comm., 2024). Larger grains were mapped with a 2-µm interval (pixel size), and smaller grains with a 1-µm interval. All Fe and P-maps obtained are presented in Figs. 4 and 5.

## Compositional profiles

To gain insights into the long-term pre-eruptive history of a volcano, the study of Fe-Mg zoning is crucial as it can record short-term pre- and syn-eruptive magmatic timescales (e.g. Bouvet de Maisonneuve et al., 2016). Diffusion chronometry modelling requires forsterite (here calculated as $Fo = Mg/(Mg + Fe + Mn + Ca + Ni)$) variations (i.e. zonation) across olivine crystals. Due to inherent crystal variability and orientation, multiple profiles perpendicular to different crystal faces within a single olivine are necessary to ensure anisotropic diffusion and statistically representative timescales (Shea et al., 2015; Appendix Fig. E1). However, acquiring hundreds of EPMA Fo profiles per sample is costly, time-consuming, and has a minimum resolution of 2 µm. An alternative is using greyscale transects on BSE images, which offer higher spatial resolution, with a minimum step size of 274 nm between points in this study. Indeed, BSE images map compositional variations based on mean atomic weight (Reed, 2005) with a strong negative linear relationship between greyscale values and Fo (Fig. 6a, b; Martin et al., 2008; Pankhurst et al., 2018; Ruth et al., 2018), where brighter pixels indicate a higher Fe concentration, and darker pixels high Mg concentrations. However, the BSE signal can also be affected by olivine's crystallographic orientation (Joy, 1974; Lloyd et al., 1987; Prior et al., 1996, 1999; e.g. see subgrain boundary in grain 1-B and 4-K, Fig. E1), necessitating a greyscale/Fo relationship definition for each grain.

First, high-resolution BSE images were obtained with the UBC EPMA, at the conditions described in the olivine section of Appendix D. EPMA BSE images were preferred over SEM images, as EPMA enables rapid acquisition of high-resolution images thanks to its high and stable current and better-quality BSE detector. Profiles were drawn on each BSE image using Adobe Illustrator (Fig. E1) to enable more efficient and precise analyses and keep a record of their locations. Greyscale profiles were then taken across each grain's BSE image with the free software ImageJ (Schneider et al., 2012) (as done by others, e.g. Hartley et al., 2016; Martin et al., 2008; Metcalfe et al., 2021; Morgan et al., 2004; Pankhurst et al., 2018; Petrone et al., 2016; Rae et al., 2016), setting the line tool to 30-pixels thick to reduce the noise and uncertainty (Morgan et al., 2004; Pankhurst et al., 2018).

For the 83 olivine grains, we analysed 246 greyscale profiles taken normal to the apparent zonation (Figs. 7, E1). Suitable transects were then carefully chosen based on best fits and following recommendations in Shea et al. (2015); i.e. avoiding merging diffusion fronts, core loss and dipping plateaux at the profile location. After removing unsuitable profiles, an average of ~3 profiles per grain (234 profiles within 81 crystals; 2 crystals are too broken to display a clear zoning) were modelled. In 77% of the grains (62/81), the profiles were perpendicular to at least two different
<figure><img src="image_4.png" /><figcaption>**Fig. 4.** X-ray maps of Group 1 olivine phenocrysts: a) Iron (Fe) maps; b) Phosphorus (P) maps (the high-P content of the melt was removed and replaced by white for clarity). c) Intepretation of the P maps. While primary and secondary branches are the first to form during olivine growth and, therefore, define its 'skeleton', dissolution fronts (dashed lines) highlight where branches have been resorbed and replaced by low-P olivine. Prim. br. is Primary branch.</figcaption></figure>

non-parallel faces, enabling verification of grain orientation and timescale consistency for each crystal (Shea et al., 2015; see Figs. 7 and Appendix E2 for timescale variability within each grain).

EPMA measurements were then conducted on at least one greyscale profile per grain to i) correlate greyscale and Fo values (Fig. 6a, b), ii) extrapolate the originally EPMA-derived diffusion transects (Fig. 6c), and iii) convert the remaining greyscale transects into Fo transects, all through the *Greyscale-to-Fo converter* code (Aufrère et al., 2024; inspired from Pankhurst et al., 2018). When multiple greyscale- and Fo-transects were obtained within one grain, the one with the best relationship ($R^2 > 0.9$) was used to convert all remaining greyscale profiles. Of the 234 greyscale profiles, 129 were chemically measured (olivine analytical conditions described in Appendix D). Profiles start in the groundmass touching the mineral rims or ~2 µm from the crystal edge, varying from 40 to 200 µm in length and with spot intervals ranging from 2 to ~4 µm. The EPMA-derived extrapolated Fo transects and the greyscale-derived-only transects were then used for diffusion chronometry (see *Diffusion Modelling* section).

# DATA TREATMENT AND MODELLING

## Pre-eruptive conditions

The rates at which Fe and Mg diffuse within olivine (and, therefore, the calculated timescales) is strongly temperature-sensitive, with a lesser influence from pressure and oxygen fugacity (Dohmen et al., 2007; Dohmen & Chakraborty, 2007a, 2007b). Therefore, these values are crucial to subsequent diffusion chronometry analysis of zoning patterns in the olivine grains.

To calculate pre-eruptive conditions, all our calculations are based on the matrix glass as its composition is the closest to that of the true pre-eruptive melt. Both the thermometer and the $fO_2$ equation used in this study require knowledge of Fe partitioning between ferrous and ferric iron; however, the glass analyses report Fe-contents in terms of $FeO_T$. Still, the parental magma Fe-partition can be constrained from the whole-rock analyses (assuming they are of the same composition), if the sample is fresh and quenched, i.e. its initial ferrous/ferric content was preserved and unchanged by later oxidation. The whole-rock analysis obtained by Wilson & Russell (2017) for sample
<figure><img src="image_5.png" /><figcaption>**Fig. 5.** X-ray maps of Group 2a olivine microphenocrysts: a) Iron (Fe) maps; b) Phosphorus (P) maps (the high-P content of the melt was removed and replaced by white for clarity). c) Interpretation of the P maps. Crystal textures range from spongy to skeletal, and show primary and some secondary branches, which are the first to form during olivine growth. All feature embayments (noted as ε) within their crystal framework, suggesting unfinished crystallisation rather than a dissolution event. P. b. is Primary branch.</figcaption></figure>

AW-15-070 returned total iron reported as FeO. The sample was thus reanalysed by ALS Global in North Vancouver (BC, Canada), with the FeO content measured volumetrically by titration (see Results). By using the whole-rock-derived Fe-partition for the glass, we i) consider the whole-rock composition as a proxy for the parental magma (i.e. assume no mixing or assimilation was involved, although a simplification), and ii) assume the Fe ratio did

not change from that of the parental magma to the pre-eruptive one (see Appendix F).

Accurate temperature assessments are key for diffusion chronometry as the appropriate diffusion coefficient (the system being treated as isothermal) is highly temperature-dependant (Dohmen & Chakraborty, 2007a, 2007b). Therefore, Putirka et al. (2007)'s equation 4 (the most appropriate olivine-liquid
<figure><img src="image_6.png" /><figcaption>Fig. 6. Example of a greyscale transect calibrated against EPMA Fo measurements (here, for gr4-A) using the Greyscale-to-Fo converter (Aufrère et al., 2024), demonstrating the Fo-greyscale negative correlation for each grain. a) Greyscale and Fo transects vs distance from the rim, b) Linear regression equation linking the filtered grey values to the Fo values, c) The Fo profile observed in a) now extended based on grey values in a) and the greyscale-to-Fo equation obtained in b).</figcaption></figure>

thermometer for hydrous melts, see Appendix G; SEE ±29°C in hydrous conditions) was used here and implemented via the open-source Python3 package Thermobar v.1.0.42 (Wieser et al., 2022). First, olivine-liquid equilibrium was tested (KD calculations, Table G2; Rhodes Diagram, Fig. 8) using in-contact

glass and olivine rim compositions, for both the phenocryst and microphenocryst populations. Pairs where KD did not fit the 0.27 to 0.33 interval (Roeder & Emslie, 1970) were then removed to enable better constraints (and smaller error) on the pre-eruptive magma temperature. Wilson & Russell (2017) determined the pre-eruptive storage conditions to be ~100 MPa, with an H2O content of 0.5 to 1 wt %. Both water-content end-members (i.e. 0.5 and 1 wt %) were considered in the temperature calculations to assess the temperature range and its consequences on timescales (maximising and minimising timescales, respectively; see Results).

The diffusion chronometry model requires the input of oxygen fugacity. To calculate it as a function of the ferric/ferrous ratio, temperature and glass composition, equation 4 of Borisov et al. (2018) was rearranged (Eq. 1) where T is in K and Xᵢ are mole fractions of the oxides. Both ferric/ferrous ratio and oxygen fugacity are expressed as log₁₀ units.

$$ (1) \log(f_{O_2}) = \frac{1}{0.207} \cdot \left[ \log\left(\frac{X_{FeO1.5}}{X_{FeO}}\right) - \frac{4633.3}{T} + 0.445 \cdot X_{SiO_2} \right. \\ \qquad \left. + 0.900 \cdot X_{TiO_2} - 1.532 \cdot X_{MgO} - 0.314 \cdot X_{CaO} \right. \\ \qquad \left. - 2.030 \cdot X_{Na_2O} - 3.355 \cdot X_{K_2O} + 4.851 \cdot X_{P_2O_5} \right. \\ \qquad \left. + 3.081 \cdot X_{SiO_2} \cdot X_{Al_2O_3} + 4.370 \cdot X_{SiO_2} \cdot X_{MgO} + 1.852 \right] $$

To report the fO₂ estimates (Eq. 1) as ΔNNO and ΔFMQ, the equations from Huebner & Sato (1970), Frost (1991) and Iacovino (2022) at 100 MPa pressure (Wilson & Russell, 2017) were used, with the temperature as constrained previously.

## Diffusion modelling

To efficiently infer timescales from Fe-Mg interdiffusion in olivine, the Excel-based AUTODIFF model was applied (e.g. Hartley et al., 2016; Pankhurst et al., 2018; Couperthwaite et al., 2020; Kahl et al., 2022, 2023a). There are two setups for AUTODIFF which correspond to either a crystal that diffuses and re-equilibrates homogeneously (continuously) from the edge, or a step function (initial square wave) between a homogeneous core and a homogeneous rim of a different composition. Here, the homogeneous crystal setup was used (Fig. 7) as profile shapes matched these models better than with the square function setup. The reader is referred to Couperthwaite et al. (2020) for a thorough description of AUTODIFF and a demo version of it, and Couperthwaite et al. (2021) and Kahl et al. (2023a) for recent upgrades to allow non-isothermal modelling of both diffusion and growth combined. Simple diffusion modelling was carried out here as the combined diffusion and growth modelling requires inputs (e.g. parental melt composition, growth, and cooling rates) that are not yet constrained due to the scarcity of geochemical studies within the MMVC.

Rim-to-core profile inputs into AUTODIFF were from either the EPMA-based-greyscale-extended profiles, or the greyscale-derived-only profiles (converted to Fo values). The initial boundary conditions that are held constant for all the grains, and for the duration of their pre-eruptive residence time are the pre-eruptive conditions (temperature, pressure, oxygen fugacity); anisotropy, FOCore and FORim are grain dependent. While the FOCore end-member composition corresponds to the beginning of the plateau, the FORim end-member (end of the steeply declining Fo value, as it is normally zoned) needed to be adjusted until the model fit was optimised. To this end, adjustments were necessary for the
<figure><img src="image_7.png" /><figcaption>**Fig. 7.** Fo contents along multiple transects oriented perpendicular to exterior crystal faces, from which maximum timescales have been obtained using the AUTODIFF model at 1097°C, 1 kbar, and NNO + 0.5. Grain gr3-B is an example illustrating how multiple traverses on a single grain give very similar chemical profiles that are well fit by AUTODIFF, along with consistent timescales within the grain that fall well within uncertainties (timescales in subscript are the lower bound of the uncertainties, and superscript are the upper bound). The quality of the traces obtained from EPMA and greyscale are equally good, and quenched rims (not considered in the diffusion model) are only apparent through the greyscale-based profiles as those extend up to the rim edge. Point symbols are larger than their related uncertainties.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 8.** Rhodes diagram (Dungan et al., 1978; Rhodes et al., 1979) showing olivine-liquid equilibrium and disequilibrium for Group 1 and 2a populations. The black solid and dotted lines represent the olivine-liquid equilibrium field defined by Roeder & Emslie (1970) ($K_D(\text{Fe-Mg})_{\text{ol-liq}} = 0.30 \pm 0.03$). Group 2a core compositions are in near-equilibrium with the averaged glass throughout the thin section. Groups 1 and 2a olivine rims are in equilibrium with the touching glass. The whole-rock composition is too mafic to be in equilibrium with Group 1 cores; removing 2 wt % of Group 1 cores from the whole rock leads to a liquid with Mg# in equilibrium with Group 1 cores (ellipse). The olivine accumulation trend is that from Fig. 3 in Putirka (2008).</figcaption></figure>

profiles displaying overgrown quenched rims (Fig. 7), as these were not considered for diffusion purposes (Couperthwaite et al., 2020).

AUTODIFF uses the corrected version (Dohmen & Chakraborty, 2007b) of the Fe-Mg diffusion coefficients obtained by Dohmen et al. (2007) and Dohmen & Chakraborty (2007a, 2007b), which accounts for composition, pressure (P), temperature (T), oxygen

fugacity ($fO_2$) and anisotropy. The latter is accounted for through crystal orientations obtained via EBSD; the P, T and $fO_2$ conditions applied here are those constrained by Wilson & Russell (2017) ($P \sim 100$ MPa; $H_2O = 0.5-1$ wt %) and in this study (see below). Minimum and maximum timescales (see Excel sheets E3, E4 for input parameters and obtained timescales) were inferred as follows: i) minimum timescales at the crystallisation temperature calculated when glass $H_2O = 0.5$ wt %, using the shortest timescale obtained within each grain (for grains in which several transects were obtained); and ii) maximum timescales, at the temperature inferred when $H_2O = 1$ wt %, using the longest timescale obtained within each grain.

# RESULTS
## Petrography

The AW-15-070 pillow rinds are vesicular (18.1 area%), glassy (sideromelane glass, 55.5 area%) and porphyritic, with 3.6 area% (5.1 wt %; see Excel sheet A2) $Fo_{33}$ olivine glomerocrysts (green in Fig. A1; Fig. 9), and 18 area% plagioclase phenocrysts (blue in Fig. A1). These olivine phenocrysts often host ~30 µm inclusions of sulphides and $Cr#_{31}$ - $Mg#_{56}$ spinel (on average, with $Cr# = 100*Cr/(Cr + Al + Fe^{3+})$ and $Mg# = 100*Mg/(Mg + Fe^{2+})$; Fig. 9 and 10). The plagioclase phenocrysts consist of large, 1-to-4-mm long sieve-textured grains (~$An_{35}$ embayed core, 6.5 area%; ~$An_{60}$ recrystallised rim) and 0.2-to-1.2 mm long, ~$An_{60}$ phenocrysts (the ~$An_{60}$ totalling 11.5 area%), with $An = 100*Ca/(Ca + Na + K + Mg)$. One highly resorbed 2-mm diameter quartz grain (dark red in Fig. A1) is also present (0.2 area%), surrounded by an augite corona (orange in Fig. A1). Microphenocrysts of 150-to-750 µm long ~$Mg#_{79}$ resorbed augite (1.1 area%) and 50-to-200 µm long $Fo_{76-79}$ olivine (3.4 area%; Fig. 9) are present in the glassy matrix, along with glomerocrysts of augite, anorthite and 60-to-400 µm long, $Fo_{75-77}$ olivine (Fig. 9). Towards the core of this pillow-lava, clusters of feathery (dendritic), quenched augite crystals mainly nucleating on plagioclase and growing outward, gradually replace the glass.
<table><thead><tr><th></th><th>Group 1</th><th>Group 2a</th><th>Group 2b</th></tr></thead><tbody><tr><td><strong>Observation</strong></td><td></td><td></td><td></td></tr><tr><td>BSE images</td><td></td><td></td><td></td></tr><tr><td>n</td><td>81 on grain-mount</td><td>20 on thin section</td><td>Studied by W &amp; R (2017)</td></tr><tr><td>L</td><td>200 µm to 1.4 mm</td><td>50 to 200 µm</td><td>60 to 400 µm</td></tr><tr><td>Classification</td><td>Phenocrysts</td><td>Microphenocrysts</td><td></td></tr><tr><td>Glomerocryst type</td><td>Olivine-only</td><td>Olivine-only</td><td>Olivine-anorthite-augite</td></tr><tr><td>Proportion</td><td>5.1 wt.%</td><td>Total: 4.8 wt.%</td><td></td></tr><tr><td>Texture</td><td>Subhedral</td><td>Subhedral to euhedral</td><td>Anhedral</td></tr><tr><td>Core-to-rim zonations</td><td>Reverse (in 1/3 of the grains) followed by normal zoning</td><td>Weak normal zoning on the most mafic grains</td><td>Unzoned</td></tr><tr><td>Core composition</td><td>Centered on Fo83</td><td>Fo76-79</td><td>Fo75-77</td></tr><tr><td>Rim composition</td><td>Fo76-80</td><td>Fo75-77</td><td>N/A</td></tr><tr><td>Inclusions</td><td>~30 µm Cr-spinel and sulphide</td><td>~5 µm spinel</td><td>Rare ~10 µm spinel</td></tr><tr><td><strong>Interpretation</strong></td><td></td><td></td><td></td></tr><tr><td>Genetic terminology</td><td>Philocrystic core, autocrystic outermost rim</td><td>Autocrysts</td><td>Autocrysts</td></tr></tbody></table>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Textural and chemical comparison of the two main olivine groups observed in the Lillooet Glacier sample, including the number of grains studied in each group (n), the length range of individual grains (L), the size classification, the glomerocryst type (which phases compose the glomerocrysts), the proportion, texture, core and rim compositions, inclusions, and genetic terminology (based on Zellmer, 2021; see Discussion section).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. 10.** Backscattered electron images of Group 1 olivine phenocrysts, featuring their dissolution features: cusps of various wavelengths and depths (a-c) as well as deep embayments (noted as ε; d-f).</figcaption></figure>
# Olivine

## Textural variations

Two olivine populations were distinguished based on size, composition, and texture: Group 1 phenocrysts, and Group 2 microphenocrysts (Fig. 9).

Group 1 olivine phenocrysts range in size from 200 µm to 1.4 mm (on average 550 µm long). These mostly occur as olivine-only glomerocrysts (aggregates composed of 2 to 20 individual grains, single crystals are uncommon; see Fig. A1) up to 2-mm long, free of deformation when observed under a petrographic microscope (i.e. they do not show undulatory extinction nor bounding kink-bands). They are subhedral, displaying irregular corners and featuring, in 36/83 grains, deep embayments or large, irregularly shaped, glassy apparent melt inclusions (Fig. 10d-f) with the same composition as the surrounding groundmass glass (Aufrère & Williams-Jones, 2024). These embayments and apparent melt inclusions intersect the inner structure (P-rich branches) of the olivine crystals, even though some have now been partially filled by low-P olivine (Fig. 4b,c). In addition, half the grains (44/83) display cuspate margins (regularly spaced hummocks and valleys; Thornber & Huebner, 1985) of variable wavelengths and depths (Fig. 10a-c). A quarter (21/83) of the grains present both cusps and embayments in their studied section (Appendix, Excel sheet H1). The grains were classified based on their apparent degree of resorption (in 2D) and whether they display embayments and/or cuspate margins, as well as how deep these are (Appendix, Excel sheet H1). Of the 83 grains, 22 were not resorbed, 23 were slightly resorbed, 21 moderately resorbed and 17 highly resorbed. No apparent correlation is seen between the level of olivine resorption and core composition (Fig. H2).

Group 2 olivine microphenocrysts are smaller (from 50 to 200-µm long) and while most (Group 2a) are subhedral to euhedral and occur both as individual grains and olivine-only glomerocrysts (Figs. A1, 5, 9), a subset (Group 2b) are rounded, not embayed, and always form part of augite-anorthite-olivine glomerocrysts (Figs. A1, 9). The latter were interpreted by Wilson & Russell (2017) to be cognate, gabbroic inclusions in the sampled lavas, that likely formed by sidewall crystallisation under pre-eruptive conditions.

Group 2a olivine microphenocrysts feature perfect angle corners (Fig. 5c); however, they also display numerous embayments and apparent melt inclusions (Fig. 5c). Indeed, half of Group 2a olivine population show embayed grains; while the most embayed grains have a spongy texture (e.g. fourth row in Fig. 5), the least embayed grains are almost euhedral (e.g. second row in Fig. 5). Phosphorus maps of those grains reveal in-or-out-of-section primary branches alongside well-defined secondary branches and embayments (in this section, appearing as large melt inclusions) centred within their structure (Fig. 5b,c). When the olivine grain is sufficiently euhedral to display sharp corners, and embayments are cross-cut by the section plane (i.e. implying that any liquid inside the olivine was in direct contact with the external melt), embayments occur between two sharp corners-apices (Fig. 5). For the spongy olivine grains, only the primary branches are observed, whereas near-euhedral grains start showing the beginning of the secondary branch framework (Fig. 5b, c). The curvature of the olivine buds and the very few crystal edges in the spongy grains suggest a 'transitional morphology between diffusion-controlled growth, which amplifies protuberances, and interface-controlled growth, which flattens faces' (Welsch et al., 2014).

## Compositional variations

The olivine populations are also chemically distinct. Group 1 olivine cores are centred on a Fo₈₃ composition and have high NiO

content (0.15–0.28 wt %; Fig. 3). Although they are not in equilibrium with the whole-rock composition (the latter being too mafic, Mg# = 61), they are in equilibrium with a whole-rock composition from which we have removed, incrementally, 2 wt % of Group 1 olivine core compositions (Table F1; Mg# = 58; Fig. 8). Group 1 olivine glomerocrysts show a continuous normal zoning (broad diffusion width up to 70 µm into the crystal) on their outer rim, which follows the embayed, cuspate margins (Figs. 4a, 10). Their rims are Fo₇₆–₈₀, display a lower NiO content (centred on 0.10 wt %; Fig. 3) and are in equilibrium with the glass (Fig. 8; Table G2). Of the 81 diffusion-suitable olivine grains, 2 distinct zoning families are noted. From the core to rim, 54 grains are normally (N) zoned, while 27 show a reverse zoning occurring before their normal zoning (RN) (Fig. 11a). Based on whether the transects were greyscale- or EPMA-based, some additional growth (G), i.e. quenched rim, can be seen on the outermost rim of most grains, for both zoning families (Figs. 7, 11a). The normal zoning of the N(G) and RN(G) families falls within the Fo₇₅–₈₄ interval. For the RN(G) family, the reverse zoning that precedes the main normal zoning is overall minor and of much smaller gradient than the normal zoning, usually spanning up to 1 Fo unit (Fig. 11), within the Fo₈₂–₈₄ interval. For this reason, it could not be modelled for timescales.

The largest and most Mg-rich Group 2a olivine (up to Fo₇₈.₇; Figs. 3, 5a) show a weak normal zoning (shorter diffusion width, up to 10 µm), while the most Fe-rich ones (up to Fo₇₅.₆; Figs. 3, 5a) are unzoned and present a spongy texture. Group 2b olivine within olivine-plagioclase-augite glomerocrysts are unzoned. Group 1 olivine phenocryst rim compositions overlap those of Group 2a microphenocryst cores, while Group 2a rims and Group 2b olivine (from the gabbroic inclusions) are of a slightly more evolved composition than Group 2a cores (up to Fo₇₅.₅; Fig. 3). Like Group 1 olivine rims, Group 2a and 2b are in equilibrium with the glass (Fig. 8; Table G2).

# Diffusion chronometry

## Pre-eruptive conditions

The whole rock FeO content obtained volumetrically is as follows: FeO = 7.88 wt %; hence, Fe₂O₃ = 1.99 wt %. These contents yield $\frac{Fe^{3+}}{\Sigma Fe} = 0.185 \pm 0.020$, which is on the lower end of the $\frac{Fe^{3+}}{\Sigma Fe}$ range for subduction zones (0.18–0.32; Kelley & Cottrell, 2009), and lower than the estimate for the nearby Mosaic Ridge (Fig. 1b) basalts (0.26–0.28; Venugopal et al., 2020). This ratio was then input into Thermobar v.1.0.42 (Wieser et al., 2022) to test for olivine-glass equilibrium for the two olivine populations. The pairs with $K_D$ that fit the 0.27–0.33 interval (Roeder & Emslie, 1970) were then used for olivine-liquid thermometry calculations, using Putirka et al. (2007)'s equation 4 assuming P = 100 MPa and H₂O = 0.5–1 wt % (Wilson & Russell, 2017; see earlier discussion). For H₂O = 0.5 wt %, Group 1 rims would have crystallised at 1106 ± 30°C and Group 2a microphenocryst rims at 1095 ± 30°C (1σ, Table G2). At slightly higher water content (H₂O = 1 wt %), the crystallisation temperature of Group 1 rims yields 1097 ± 30°C, and that of Group 2a rims 1086 ± 30°C (1σ; Table G2). For diffusion chronometry purposes, the important temperature is that of the Group 1 olivine rims; therefore, both 1106 ± 30°C and 1097 ± 30°C were used as the temperature at which diffusion occurred. Using Eq. 1 (equation 4 of Borisov et al., 2018), both water-dependant, pre-eruptive temperatures yield the same oxygen fugacity conditions: $\Delta NNO = 0.5 \pm 1.1$ ($\Delta FMQ = 1.2 \pm 1.1$). This is comparable to the NNO + 0.74 conditions inferred from olivine-hosted melt inclusions from Mosaic Ridge (Fig. 1b; Venugopal et al., 2020). Additionally, it lies within the $\Delta FMQ$ range (0.6–1.3) constrained at Cracked Mountain, just south of the Mount Meager massif (Fig. 1b; Harris & Russell, 2022).
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** a) Examples of representative Fo profiles (shown as distance from rim) for each Group 1 zoning family: Normal (—Growth), N(G), and Reverse-Normal (—Growth), RN(G). b) Violin plots showing the Fo-content of the core (profile plateau) for both the N(G) and RN(G) families, and the maximum Fo obtained at the peak of the reverse zoning, for the RN(G) family. Note the lower Fo core content of the RN(G) family compared to the N(G) one, and the overlap between the Fo core of the N(G) family and the Fo max of the RN(G) family.</figcaption></figure>

## Timescales

Using AUTODIFF and the pre-eruptive conditions mentioned above (100 MPa, 1106–1097°C, NNO + 0.5), 234 timescales (for each temperature) were calculated for a total of 81 diffusion-suited Group 1 olivine grains, based on their primary, normal zoning (Figs. 7, 11a, E1) throughout both the N(G) and RN(G) populations.

The 81 representative timescales (one per grain) were then compared based on their zoning families; each calculated timescale has an average uncertainty of 0.30 log units (1σ). Group 1 N(G) and RN(G) yielded very similar timescales with a minimum peak (for the Kernel density estimation plots) at 27 and 23 days, respectively (Fig. 12a), and a maximum peak at 56 and 74 days (Fig. 12b). Thus, considered together, 68% of the minimum timescales plot within 1.5 months (with a peak at 26 days; Fig. 12a), and 68% of the maximum timescales within 3.5 months (with a peak at 56 days; Fig. 12b), suggesting a very restricted period of residence in the magma reservoir under pre-eruptive conditions.

## DISCUSSION

Our investigation into olivine populations within the LG basalts yields crucial insights into the magmatic processes preceding the eruption. Firstly, distinct olivine groups, characterised by their texture and composition, are discerned. Secondly, thermometry and fO₂ calculations enable further constraints on pre-eruptive conditions, building on Wilson & Russell (2017)'s pressure and water-content estimates. Lastly, Fe-Mg interdiffusion-based timescales reveal short-lived magmatic processes before eruption, with most falling within 1 to 3 months. However, before interpreting these timescales, it is important to evaluate what caused Group 1 olivine to be zoned. Interestingly, the observed normal zoning in Group 1 olivine is parallel to their embayments and cuspat margins (Figs. 4a, 10), indicating that these occurred before the olivine became normally zoned. Understanding the magmatic processes that preceded the LG eruption and caused such textural and compositional variations in the olivine is important, as it contributes to eruption forecasting initiatives within the MMVC.

### Sampling and transport of Fo₈₃ olivine crystal mush

Most of Group 1 olivine phenocrysts occur as glomerocrysts that are not tightly interlocked but rather tend to show pore space between crystal faces (Fig. A3a, b). Aside from their embayed

crystal faces, they display irregular corners and oxides trapped along the edges (e.g., Fig. 10b). These olivine crystals are also free of deformation features (such as the bounding kink-bands observed in Helz's, 1983, 1987 Class 1 grains, interpreted as cumulates) when observed under a petrographic microscope. These observations are consistent with a magmatic accumulation of crystals, i.e. phenocrysts forming a 'crystal mush' (e.g. Marsh, 1996; Jerram et al., 2003) through, for example, sidewall crystallisation. Moreover, Group 1 olivine glomerocrysts and isolated crystals share similar textures and composition, so are interpreted here as being fragments of the remobilised mush (Fig. 13; Burgisser & Bergantz, 2011; Hughes et al., 2021). The whole-rock composition is also too Mg-rich to have crystallised Group 1 olivine (Mg# = 61; Fig. 8) and results from olivine accumulation. However, if the whole-rock composition is adjusted by subtraction of only ~2 wt % olivine (i.e. Group 1 olivine core compositions), the new melt composition (Mg# = 58; Fig. 8) could have been in equilibrium with Group 1 olivine cores. While there is a higher (~5 wt %) content of Group 1 olivine in thin section (Fig. A1), the thin section or the whole rock composition may not be representative of the system (i.e. perhaps anomalously enriched or depleted in olivine). Nevertheless, 2 to 5 wt % is comparable at this scale and supports olivine accumulation. These features support the hypothesis whereby Group 1 olivine originate from a reworked crystal mush and were picked up during ascent, ultimately increasing the whole rock Mg#. Geophysical studies reveal a potential deep magmatic system at ~5–15 km below sea level beneath the Mount Meager Volcanic Complex (Hanneson & Unsworth, 2023), and one petrological study conducted on Mosaic Ridge (Venugopal et al., 2020) suggests polybaric olivine crystallisation ranging from 100 to 400 MPa, i.e. ~3–12 km. Therefore, the reworked crystal mush may have originated from a similar depth range of ~3–15 km (Fig. 13).

The apparent melt inclusions are large, irregularly shaped, and glassy (Fig. 10d-f) as opposed to being ellipsoidal or faceted and equant, but also have the same composition as the surrounding glass (Aufrère & Williams-Jones, 2024). Therefore, they more likely represent embayments (connected to the surrounded glass 'in front' or 'behind' the olivine cross-section) than primitive melt inclusions that sampled olivine's parental magma during their growth. About half the Group 1 olivine grains (36/83) are embayed (Fig. 10d-f; Excel Sheet H1); however, embayments can either be a sign of a dissolution event or the result of rapid growth (Welsch et al., 2013). To examine the cause of
<figure><img src="image_14.png" /><figcaption>**Fig. 12.** a) Minimum and b) Maximum timescales inferred from Fe-Mg diffusion chronometry in 81 olivine crystals, plotted against their crystal index i.e. their rank when sorting timescales from smallest to largest (e.g. Hartley *et al.*, 2016; Pankhurst *et al.*, 2018; Rasmussen *et al.*, 2018). 68% of the calculated timescales fall within 1.5 and 3.5 months. Insets are Kernel density estimation plots reflecting the density of the minimum and maximum timescales obtained for each of the zoning families. All grains return similar timescale peaks (even between the two zoning families), so that the overall trend, N(G) + RN(G), shows a minimum peak at 26 days and a maximum peak at 56 days. Although all form a continuum, 54 grains belong to the N(G) population and 27 to the RN(G) one, hence the N(G) + RN(G) trend being weighted towards the N(G) population.</figcaption></figure>

embayments in Group 1 olivine, P-maps were made of these phenocrysts and revealed P-rich branches (i.e. the olivine 'skeleton') cross-cut by low-P olivine (Fig. 4b,c). This indicates that these branches were partially resorbed by a melt pool and gradually replaced by newly-crystallising low-P olivine, although they did not have time to fill the voids completely (as seen

in Bouvet de Maisonneuve *et al.*, 2016). Interestingly, these P-maps do not reveal several resorption events (such as that observed at Kilauea; Mourey *et al.*, 2022), but rather only one.

Moreover, the cuspate margins observed in more than half of the Group 1 olivine grains (44/83; Fig. 10a-c) are similar to those produced experimentally by Thornber & Huebner (1985; their
<figure><img src="image_15.png" /><figcaption>**Fig. 13.** Schematic cross-section of the Lillooet Glacier transcrustal magmatic system based on olivine phenocrysts and microphenocrysts observed in the Lillooet Glacier pillow-lavas (Group 1, 2a and 2b; Fig. 9). Hot, mafic magma inputs at the base of the magmatic system likely entrained Group 1 olivine phenocrysts on their way up, causing resorption and for some, development of reverse zoning, prior to settling in the pre-eruptive magma chamber. While Group 2b olivine crystallised on the magma chamber margins, Group 1 developed rims of Fo₇₆ composition, and skeletal Group 2a microphenocrysts started to crystallise. During the following approximately 1 to 3 months, diffusion occurred to re-homogenise the core-rim boundary of Group 1 olivine, until it stopped upon eruption.</figcaption></figure>

Fig. 8a-c, f-h) while studying dissolution of olivine in basaltic liquids, further supporting a resorption event. For olivine grains incorporated in a superheated (relative to the entrained olivine) basaltic melt, Thornber & Huebner (1985) inferred that the wavelength of cusps increases with decreasing temperature of the melt. On the other hand, when resorption occurs in a silica-enriched liquid, the olivine crystals are first rounded, and with decreasing temperature, shorter wavelength cusps form. Here, Group 1 olivine are not rounded but rather display cuspatе margins of varying wavelength (from 10 to 160 µm; Fig. 10a-c). Therefore, prior to reaching the pre-eruptive reservoir, they were most likely incorporated in a hot, basaltic magma and as this melt cooled and approached the olivine-saturation surface, the resorption rate slowed with a resulting increase in cusp wavelength (Thornber & Huebner, 1985). These cuspatе margins, therefore, record a second hot, basaltic magma injection.

For the 83 olivine grains, no correlation was found between Fo content and the degree of resorption (Fig. H2). One hypothesis, based on the Helz (1987) model, would be that their degree of resorption depends on how the olivine crystals were entrained in their new host melt. In this scenario, a thin layer of liquid surrounding the olivine crystals became superheated during entrainment and was no longer in equilibrium with the crystals. This resulted in the (at least partial) resorption of the olivine grains, until the adjacent melt reached equilibrium with the olivine. If the adjacent melt of one grain continued to be supplied by

more mafic hot magma, this grain would likely undergo stronger dissolution. Such a model would explain not only the local resorption experienced by some olivine grains, but also the fact that the degree of resorption was independent of the olivine composition at the time (now its core composition). Moreover, various degrees of resorption reflect a continuous process of entrainment and continual mush disaggregation: as olivine glomerocrysts are entrained into a melt, the crystals on their outside become heavily embayed and ultimately detach, exposing new olivine to the melt for dissolution.

Reverse zoning, observed in a third of the Group 1 (the RN(G) family; Fig. 11a) olivine phenocrysts, further supports the hypothesis whereby these olivine grains were entrained in a hot, basaltic melt. Indeed, based on their relative compositions, it appears the N(G) and RN(G) Group 1 olivine crystallised from the same parental melt (Fig. 8), although the core content of the RN(G) grains is ~1 Fo unit less forsteritic than that of the N(G) family (Fig. 11b). Such reverse zoning developed post-dissolution, immediately after equilibrium was attained, so that the exterior of the crystals re-equilibrated to the new melt Fo equilibrium. The reason why not all grains recorded such reverse zoning could be that the RN(G) grains recorded the earliest, hotter stages of the magma mixing, and/or had the lowest Fo cores. Indeed, the reverse zoning Fo overlaps the plateau of the N(G) family (~Fo₈₃; Fig. 11b), so a mafic input in equilibrium with Fo₈₃ would have not caused a reverse zoning in the already-Fo₈₃ crystals. Still, it is
worth noting that the mafic input causing this reverse zoning was of the same Mg# as the parental magma of Group 1 olivine (to be in equilibrium with Fo₃₃ crystals in these conditions, liquid Mg# = 58; Fig. 8). Therefore, at least two Mg#₅₈ magma inputs occurred: one from which the Group 1 olivine phenocrysts crystallised, and another causing a third of them to be reversely zoned. Following this scenario, the N(G) and RN(G) Group 1 olivine would have been entrained in a hot, mafic magma (Fig. 13).

As Zellmer (2021)'s definitions apply to part of or to a whole crystal, Group 1 olivine cores can thus be considered philocrystic (antecrysts related to a broad transcrustal magmatic system; Mangler et al., 2020; Figs. 9, 13). Such processes, in which philocrysts are entrained and partially 'digested' by a recharge magma on its way to the surface has been dubbed 'petrological cannibalism' (Cashman & Blundy, 2013). Similar processes of magma recharge, crystal mush remobilisation and mixing within a polybaric, transcrustal magmatic system have been described in many volcanic arcs, including the Cascades: Mount Rainier (Stockstill et al., 2003), Mount St. Helens (Leeman & Smith, 2018) and Mount Hood (Koleszar et al., 2012).

## Entrainment of gabbroic and granitic minerals

Gabbroic inclusions as well as quartz and sieve-textured, albite-rich plagioclase phenocrysts are observed in the LG basalts; these were likely dislodged and entrained either by the rising pulses of hot mafic magma or upon eruption, during magma ascent. Indeed, the unzoned Group 2b olivine microphenocrysts only occur as part of olivine-anorthite-augite glomerocrysts and are of slightly more evolved composition than Group 2a olivine (down to Fo₇₅; Figs. 3, 9). Olivine crystallised in the cooler parts of a magma reservoir may display a higher fayalitic composition compared to those in equilibrium with the melt. This observation aligns with the sidewall crystallisation theory proposed by Wilson & Russell (2017) (Fig. 13).

Additionally, large, sieve-textured plagioclase phenocrysts with ~An₃₅ cores (light blue; Fig. A1) and ~An₆₀ rims (dark blue; Fig. A1), along with resorbed, augite-corona-bearing quartz phenocrysts (red; Fig. A1) were observed. The core of these plagioclase and quartz phenocrysts are not in equilibrium with their carrier melt; therefore, they likely represent recycled philocrysts that grew autocrystic reaction rims, upon their assimilation in their carrier melt. Two hypotheses for their origin are proposed: either they are the result of country rock assimilation (as basement rocks at the MMVC are of granodiorite composition), or they originate from a felsic crystal mush beneath the MMVC and were picked up during magma rise upon eruption.

## Rapid growth of the Fo₇₆ olivine population prior to eruption

Group 1 olivine rim compositions overlap those of Group 2a microphenocryst cores and rims (Fig. 3), suggesting they grew concurrently from a same melt in equilibrium with Fo₇₆. The Rhodes Diagram (Fig. 8) and K<sub>D</sub> (Table G2) show that these olivine groups are in equilibrium with the glass. Therefore, this ~Fo₇₆ population crystallised immediately prior to eruption and is autocrystic (Fig. 9).

Phosphorus maps reveal that i) Group 2a olivine are skeletal (Fig. 5b, c), hence their embayed texture and ii) P-rich branches in Group 1 rims appear to distinguish the Fo-rich core from the Fo-poor rims (Fig. 4b, c). This suggests rapid crystallisation of a more-fayalitic Group 1 olivine overgrowth to equilibrate with the new host melt (Dohmen et al., 2017), likely around the same time as Group 2a olivine crystallisation.

Textural evidence further supports this scenario. Group 2a olivine microphenocrysts show embayments on prominent crystal faces, rather than at the junctions of crystal faces (i.e. corners; Fig. 5). This is important because, if the embayments were dissolution-related, the crystal corners and edges would have been affected first and become rounded (Noiriel et al., 2020). Indeed, the Gibbs free energy is greater for a sharp corner than for a flat surface, which makes corners more soluble than any other part of the olivine in a basaltic solvent (Mourey et al., 2022). Moreover, crystal corners diffuse into a larger melt volume (a cone) than the centre of a crystal face which can only diffuse outward effectively linearly perpendicular to the face. Therefore, diffusion-mediated dissolution will corrode corners fastest. However, here, the olivine corners are sharp with perfect olivine angles (e.g. 81° and 140° in the b-c plane, see Fig. 5c) as highlighted in Fig. 2 of Welsch et al. (2013). Both compositional and textural evidence thus suggest rapid and incomplete crystallisation, prior to eruption.

Group 2a olivine microphenocrysts were originally thought to be the same as Group 2b crystals and representing the breakup of gabbroic inclusions (Wilson & Russell, 2017). However, after this more focused study, it is now clear that Group 2a consists of morphology-transitioning olivine with most stuck in the dendritic/skeletal stage, having not had sufficient time to mature (Welsch et al., 2014; Mourey et al., 2022). Indeed, some of the Group 2a olivine are slightly more evolved (Fo₇₅) than others in the 2a population and are characterised by unzoned, spongy olivine (Fig. 5b, c). The most mafic Group 2a olivine also developed slight normal zoning, with their rims of similar composition (Fig. 5a). Likewise, normally zoned, quenched rims (Figs. 7, 11) of the same composition are present in most Group 1 olivine (given their profile is greyscale-based) in the last 5 µm of their outermost rim. As these did not have time to chemically relax (quenched rims still show a sharp compositional boundary), they must have formed at the top of the conduit where cooling rates were at their highest, and immediately prior to eruption (Pankhurst et al., 2018), at between 1086 and 1095°C (±30°C at 1σ; Table G2).

## Time under pre-eruptive conditions

Based on the core-rim delineation observed in Group 1 P-maps (Fig. 4b, c), the Fo₇₆ rims appear to have started by a phase of rapid growth. Therefore, the diffusion observed and studied for diffusion chronometry probably started at the melt-core interface but finished at a rim-core level (Dohmen et al., 2017). Through diffusion chronometry carried out on 234 normal-zoning profiles across 81 crystals from Group 1 population, 68% of the minimum timescales fall within 1.5 months (Fig. 12a), and maximum timescales within 3.5 months (Fig. 12b). These scenarios reveal a timescale peak at 26 and 56 days, respectively (Fig. 12a, b). This 1-to-3-month period reflects Group 1 philocryst residence time in a shallow (100 MPa; ~3 km depth) and more evolved reservoir, prior to eruption, where they developed autocrystic rims (Figs. 9, 13). These timescales are comparable but shorter than the 5 months to 6.3 years of residence time in an evolving magma estimated for Cinder Cone in Lassen Volcanic National Park (southernmost centre of the Cascade Volcanic Arc; Walowski et al., 2019). Our results also overlap the mixing-to-eruption timescales obtained on the Blue Lake Crater olivine (~10–60 days; Johnson & Cashman, 2020) and the Mount Hood plagioclase (a few weeks or less; Kent et al., 2010), which are the shortest reported timescales amongst the Cascade Volcanic Arc.
## CONCLUSION & PERSPECTIVES

Thorough petrological (compositional and textural) work conducted on olivine populations within the Lillooet Glacier basalts complements previous studies (Wilson & Russell, 2017; Venu-gopal et al., 2020) and informs on the dynamics and architec-ture of crustal magma reservoirs within the Canadian Cascade Volcanic Arc. Indeed, below the LG volcanic centre, a deep (~3–15 km) and mafic magma reservoir in equilibrium with Fo₈₃ olivine (Group 1) was supplied by a new injection of hot magma of the same Mg#, which passed through an olivine crystal mush (Fig. 13). This melt likely reheated the crystal mush from below, broke apart some of the olivine aggregates (Burgisser & Bergantz, 2011) and transported them into shallower melt horizons, while also resorbing them (Figs. 10, 13). Once closer to equilibrium with this Fo₈₃ olivine composition, the most evolved olivine glome-rocrysts started developing reverse zoning of Fo₈₃ composition (Figs. 11, 13). Group 1 olivine phenocrysts appear to have experi-enced two resorption events and thus were likely remobilised at least twice, prior to eruption.

The magma entraining Group 1 olivine philocrysts ultimately reached a shallow (~3 km depth), colder and more evolved pre-eruptive magma reservoir, which promoted growth of a more evolved, Fo₇₆ population (both as rims for Group 1 olivine, and as skeletal Group 2a microphenocrysts; Figs. 8, 9), and Fo₇₅ Group 2b sidewall crystallisation (Fig. 13).

Pre-eruptive conditions (P, T, H₂O, fO₂) constrained by Wilson & Russell (2017) and in this study allow us to conduct diffusion chronometry on Group 1 population and yield pre-eruptive resi-dence timescales of only 1 to 3 months (Fig. 12a, b). Previous stud-ies of products of recent eruptions (e.g. Mt. Etna, 1991–1993 and 2006, Kîlauea, 2018 and Stromboli, 2020) have shown olivine zon-ing and diffusion timescales to be correlated with eruption trig-gers (e.g. via seismicity, surface deformation or gas ratio changes; Kahl et al., 2011, 2013; Mourey et al., 2023; Voloschina et al., 2023). Within the Canadian segment of the Cascade Volcanic Arc, detailed petrochemical studies of Pleistocene and Recent eruption products are lacking and monitoring networks are embryonic (Kelman & Wilson, 2024). As such, providing an estimate of pre-eruptive primers and timescales is a significant step forward. Moreover, the LG eruption is not the only mafic eruption at the MMVC (Fig. 1b), and while many were glaciovolcanic in nature (Lil-looet Glacier, Wilson & Russell, 2017, 2018; Eastern Lillooet Ridge and Lillooet River; Harris et al., 2023; Cracked Mountain, Harris et al., 2023; Read, 1990); future mafic eruptions in the modern deglaciated landscape would not be impounded by ice (e.g. the Cheakamus basalts; Borch et al., 2023). Thus, future eruptions of low viscos-ity magmas could travel significantly greater distances with an increased likelihood of impacting populations and infrastructure.

Subsequent petrological studies should aim to investigate additional mafic and felsic eruptions at the MMVC, including geobarometry and melt inclusion volatile analyses, to better constrain crystallisation depths, eruption primers, triggers and timescales and understand the magmatic system as a whole. Given the limited existing knowledge of its volcanic history and the nascent monitoring, constraining magma residence time provides important insights for future eruption forecasting.

## ACKNOWLEDGEMENTS

We respectfully acknowledge that the area that inspired this research is within the traditional unceded territories of the Lîlwat First Nation. We thank Alexander Wilson (University of British Columbia) for access to the UBC Lillooet Glacier samples and

Christophe Névado (Université de Montpellier) for his help prepar-ing the grain-mounts. We would like to thank Fabrice Barou (Université de Montpellier) for his assistance with the EBSD, Emmy Voyer (Laboratoire Magma et Volcans) with the SEM, and Jean-Luc Devidal (Laboratoire Magma et Volcans) and Anette von der Handt (University of British Columbia) with the EPMAs. SA thanks Yannick Le Moigne (Geological Survey of Canada) for the helpful discussions about magmatic processes. We would also like to thank Juan Anzieta (Simon Fraser University) for his coding skills that allowed the development of the *Greyscale-to-Fo* Python code. Lastly, we are thankful to our three anonymous reviewers and editor, A. Kent, for thoughtful comments that greatly improved this manuscript.

## FUNDING

This research was supported by Natural Sciences and Engineering Research Council of Canada Discovery grants (RGPIN-2016-03953, RGPIN-2023-04284) to GWJ.

## DATA AVAILABILITY

Analytical data for this study is available in the EarthChem. org database (https://doi.org/10.60520/IEDA/113148) and Python *Greyscale-to-Fo* converter code via the Zenodo.org database (https://doi.org/10.5281/zenodo.10807304).

## SUPPLEMENTARY DATA

Supplementary data are available at *Journal of Petrology* online.

## REFERENCES

* Andrews, G. D. M., Russell, J. K. & Stewart, M. L. (2014). The history and dynamics of a welded pyroclastic dam and its failure. *Bulletin of Volcanology* **76**(4), 1–16. https://doi.org/10.1007/s00445-014-0811-0.

* Aufrère, S. M. & Williams-Jones, G. (2024). Glass, olivine, spinel, clinopyroxene and plagioclase analyses from the Lillooet Glacier basalts, Mount Meager Volcanic Complex (British Columbia, Canada), version 1.0. *Interdisciplinary Earth Data Alliance (IEDA)*, 811. https://doi.org/10.60520/IEDA/113148.

* Aufrère, S. M., Anzieta, J. & Williams-Jones, G. (2024). *Greyscale-to-Fo* converter for olivine diffusion chronometry. *Zenodo*. https://doi.org/10.5281/zenodo.10807304.

* Bachmann, F., Hielscher, R. & Schaeben, H. (2010). Texture analysis with MTEX-free and open-source software toolbox. *Solid State Phe-nomena* **160**, 63–68. https://doi.org/10.4028/www.scientific.net/SSP.160.63.

* Borch, A., Russell, J. K. & Barendregt, R. (2023). Cheakamus basalt lavas, British Columbia: a Pleistocene record of rapid, continu-ous eruption within a mountainous drainage system. *Canadian Journal of Earth Sciences* **60**(11), 1544–1572. https://doi.org/10.1139/cjes-2023-0004.

* Borisov, A., Behrens, H. & Holtz, F. (2018). Ferric/ferrous ratio in silicate melts: a new model for 1 atm data with special emphasis on the effects of melt composition. *Contributions to Mineralogy and Petrology* **173**(12). https://doi.org/10.1007/s00410-018-1524-8.

* Bouvet de Maisonneuve, C., Costa, F., Huber, C., Vonlanthen, P., Bachmann, O. & Dungan, M. A. (2016). How do olivines record magmatic events? Insights from major and trace element zoning. *Contributions to Mineralogy and Petrology* **171**(6), 56. https://doi.org/10.1007/s00410-016-1264-6.
du Bray, E. A. & John, D. A. (2011). Petrologic, tectonic, and met-
allogenic evolution of the ancestral cascades magmatic arc,
Washington, Oregon, and northern California. Geosphere 7(5),
1102–1133. https://doi.org/10.1130/GES00669.1.

Burgisser, A. & Bergantz, G. W. (2011). A rapid mechanism to remo-
bilize and homogenize highly crystalline magma bodies. Nature
471(7337), 212–215. https://doi.org/10.1038/nature09799.

Cashman, K. & Blundy, J. (2013). Petrological cannibalism: the chem-
ical and textural consequences of incremental magma body
growth. *Contributions to Mineralogy and Petrology* **166**(3), 703–729.
https://doi.org/10.1007/S00410-013-0895-0.

Cashman, K. V., Sparks, R. S. J. & Blundy, J. D. (2017). Verti-
cally extensive and unstable magmatic systems: a unified view
of igneous processes. *Science* **355**(6331). https://doi.org/10.1126/
science.aag3055.

Clague, J. J. & Ward, B. (2011). Chapter 44 - Pleistocene glaciation of British Columbia. In: Ehlers, J., Gibbard P. L. & Hughes P. D. (eds) *Quaternary Glaciations – Extent and Chronology: A Closer Look, Developments in Quaternary Science*, 15. Amsterdam, The Netherlands: Elsevier, pp.563–573. https://doi.org/https://doi.org/10.1016/B978-0-444-53447-7.00044-1.

Clague, J. J., Evans, S. G., Rampton, V. N. & Woodsworth, G. J. (1995). Improved age estimates for the White River and Bridge River tephras, western Canada. *Canadian Journal of Earth Sciences* **32**(8), 1172–1179. https://doi.org/10.1139/e95-096.

Connelly, J. P., Sepúlveda, S. A. & Williams-Jones, G. (2024). Qwelqwelústen / mount meager volcanic complex, British Columbia: inter-eruptive landslide susceptibility assessment using statistical machine learning techniques. *Volcanica, Accepted*.

Couperthwaite, F. K., Thordarson, T., Morgan, D. J., Harvey, J. & Wilson, M. (2020). Diffusion timescales of magmatic processes in the Moinui lava eruption at Mauna Loa, Hawaii, as inferred from bimodal olivine populations. *Journal of Petrology* **61**(7). https://doi.org/10.1093/petrology/egaa058.

Couperthwaite, F. K., Morgan, D. J., Pankhurst, M. J., Lee, P. D. & Day, J. M. D. (2021). Reducing epistemic and model uncertainty in ionic inter-diffusion chronology: a 3D observation and dynamic modeling approach using olivine from piton de la Fournaise, La Réunion. *American Mineralogist* **106**(3), 481–494. https://doi.org/10.2138/am-2021-7296CCBY.

Dohmen, R. & Chakraborty, S. (2007a). Fe–Mg diffusion in olivine
II: point defect chemistry, change of diffusion mechanisms and
a model for calculation of diffusion coefficients in natural
olivine. *Physics and Chemistry of Minerals* **34**(6), 409–430. https://
doi.org/10.1007/S00269-007-0158-6.

Dohmen, R. & Chakraborty, S. (2007b). Fe–Mg diffusion in olivine II: point defect chemistry, change of diffusion mechanisms and a model for calculation of diffusion coefficients in natural olivine. *Erratum. Physics and Chemistry of Minerals* **34**(8), 597–598. https://doi.org/10.1007/s00269-007-0185-3.

Dohmen, R., Becker, H. W. & Chakraborty, S. (2007). Fe–Mg diffusion in olivine I: experimental determination between 700 and 1,200°C as a function of composition, crystal orientation and oxygen fugacity. *Physics and Chemistry of Minerals* **34**(6), 389–407. https://doi.org/10.1007/S00269-007-0157-7.

Dohmen, R., Faak, K. & Blundy, J. D. (2017). Chronometry and speedometry of magmatic processes using chemical diffusion in olivine, plagioclase and pyroxenes. *Reviews in Mineralogy and Geochemistry* **83**(1), 535–575. https://doi.org/10.2138/rmg.2017.83.16.

Dungan, M. A., Long, P. E. & Rhodes, J. M. (1978). Magma mixing at mid-ocean ridges: evidence from legs 45 and 46-DSDP. *Geophysical

Research Letters 5(6), 423–425. https://doi.org/10.1029/GL005I006P00423.

Easterbrook, D. J. (1975). Mount Baker eruptions. *Geology* **3**(12), 679–682. https://doi.org/10.1130/0091-7613(1975)3<679:MBE>2.0.CO;2

Frost, B. R. (1991). Introduction to oxygen fugacity and its petrologic importance. *Reviews in Mineralogy and Geochemistry* **25**(1), 1–10. https://doi.org/10.1515/9781501508684-004.

Ghiorso, M. S. & Gualda, G. A. (2015). An H₂O–CO₂ mixed fluid saturation model compatible with rhyolite-MELTS. *Contributions to Mineralogy and Petrology* **169**, 1–30. https://doi.org/10.1007/s00410-015-1141-8.

Green, N. L. (1977). Multistage andesite genesis in the Garibaldi Lake area, southwestern British Columbia [PhD Thesis. University of British Columbia, Vancouver (BC, Canada)]. http://hdl.handle.net/2429/21529

Green, N. L., Armstrong, R. L., Harakal, J. E., Souther, J. G. & Read, P. B. (1988). Eruptive history and K-Ar geochronology of the late Cenozoic Garibaldi volcanic belt, southwestern British Columbia. *Geological Society of America Bulletin* **100**(4), 563–579. https://doi.org/10.1130/0016-7606(1988)100<0563:EHAKAG>2.3.CO;2

Gualda, G. A. R., Ghiorso, M. S., Lemons, R. V. & Carley, T. L. (2012). Rhyolite-MELTS: a modified calibration of MELTS optimized for silica-rich, fluid-bearing magmatic systems. *Journal of Petrology* **53**(5), 875–890. https://doi.org/10.1093/petrology/egr080.

Hanneson, C. & Unsworth, M. J. (2023). Magnetotelluric imaging of the magmatic and geothermal systems beneath Mount Meager, southwestern Canada. *Canadian Journal of Earth Sciences* **60**(10), 1385–1403. https://doi.org/10.1139/CJES-2022-0136.

Harris, M. A. & Russell, J. K. (2022). Polymagmatic Glaciovolcanism: Cracked Mountain Tuya, Canadian Cascades. *Frontiers in Earth Science* **10**, 409. https://doi.org/10.3389/feart.2022.859794.

Harris, M. A., Muhammad, M., Williams-Jones, G. & Russell, J. K. (2020) Bedrock mapping for Mount Meager Geothermal Research Initiative. In: *Geoscience BC Report (Chapter 2)*.

Harris, M. A., Russell, J. K., Muhammad, M. & Williams-Jones, G. (2022). Mount meager volcanic complex, Garibaldi Volcanic Belt, British Columbia: expanded bedrock map including Cracked Mountain, north Lillooet ridge, and west mount meager. *Geological Survey of Canada, Open File* **8881**. https://doi.org/10.4095/329886.

Harris, M. A., Russell, J. K., Wilson, A. & Jicha, B. (2023). A 500 ka record of volcanism and paleoenvironment in the northern Garibaldi Volcanic Belt, British Columbia. *Canadian Journal of Earth Sciences* **60**(4), 401–421. https://doi.org/10.1139/CJES-2022-0101.

Hartley, M. E., Morgan, D. J., Maclennan, J., Edmonds, M. & Thor-
darson, T. (2016). Tracking timescales of short-term precursors
to large basaltic fissure eruptions through Fe–Mg diffusion in
olivine. *Earth and Planetary Science Letters* **439**, 58–70. https://doi.
org/10.1016/j.epsl.2016.01.018.

Helz, R. T. (1983). Diverse olivine populations in lavas of the 1959 eruption of Kilauea volcano, Hawaii. *EOS Transaction of the American Geophysical Union* **64**, 900.

Helz, R. T. (1987) Diverse olivine types in lava of the 1959 eruption of Kilauea Volcano and their bearing on eruption dynamics. Volcanism in Hawaii. In: Decker R. W., Wright T. L. & Stauffer P. H. (eds) U.S. Geological Survey Professional Paper, 1350, pp.691–722.

Hickson, C. J., Russell, J. K. & Stasiuk, M. V. (1999). Volcanology of the 2350 B.P. Eruption of Mount Meager Volcanic Complex, British Columbia, Canada: implications for hazards from eruptions in topographically complex terrain. *Bulletin of Volcanology* **60**(7), 489–507. https://doi.org/10.1007/s004450050247.
Hollyday, A. E., Leiter, S. H. & Walowski, K. J. (2020). Pre-eruptive storage, evolution, and ascent timescales of a high-Mg basaltic andesite in the southern Cascade Arc. Contributions to Mineralogy and Petrology **175**(9), 88. https://doi.org/10.1007/s00410-020-01730-z.

Huebner, J. S. & Sato, M. (1970). The oxygen fugacity-temperature relationships of manganese oxide and nickel oxide buffers. *American Mineralogist* **55**(5-6), 934-952.

Hughes, G. E., Petrone, C. M., Downes, H., Varley, N. R. & Hammond, S. J. (2021). Mush remobilisation and mafic recharge: a study of the crystal cargo of the 2013–17 eruption at Volcán de Colima, Mexico. *Journal of Volcanology and Geothermal Research* **416**, 107296. https://doi.org/10.1016/j.jvolgeores.2021.107296.

Iacovino, K. (2022). Calculate fO2 buffer (1.5). *Zenodo*. https://doi.org/10.5281/zenodo.7387196.

Iddings, J. P. (1892). On the crystallization of igneous rocks. *Bulletin of the Philosophical Society of Washington* **11**, 71–112.

Jerram, D. A. & Martin, V. M. (2008). Understanding crystal populations and their significance through the magma plumbing system. *Geological Society, London, Special Publications* **304**(1), 133–148. https://doi.org/10.1144/SP304.

Jerram, D. A., Cheadle, M. J. & Philpotts, A. R. (2003). Quantifying the building blocks of igneous rocks: are clustered crystal frameworks the foundation? *Journal of Petrology* **44**(11), 2033–2051. https://doi.org/10.1093/petrology/egg069.

Johnson, E. R. & Cashman, K. V. (2020). Understanding the storage conditions and fluctuating eruption style of a young monogenetic volcano: Blue Lake crater (<3 ka), High Cascades, Oregon. *Journal of Volcanology and Geothermal Research* **408**, 107103. https://doi.org/10.1016/j.jvolgeores.2020.107103.

Joy, D. C. (1974) Electron channeling patterns in the SEM. In: Holt D. B., Muir M. D., Boswarva I. M. & Grant P. R. (eds) *Quantitative Scanning Electron Microscopy*. New York: Academic Press.

Kahl, M., Chakraborty, S., Costa, F. & Pompilio, M. (2011). Dynamic plumbing system beneath volcanoes revealed by kinetic modeling, and the connection to monitoring data: an example from Mt. Etna. *Earth and Planetary Science Letters* **308**(1–2), 11–22. https://doi.org/10.1016/j.epsl.2011.05.008.

Kahl, M., Chakraborty, S., Costa, F., Pompilio, M., Liuzzo, M. & Viccaro, M. (2013). Compositionally zoned crystals and real-time degassing data reveal changes in magma transfer dynamics during the 2006 summit eruptive episodes of Mt. Etna. *Bulletin of Volcanology* **75**(2), 692. https://doi.org/10.1007/s00445-013-0692-7.

Kahl, M., Morgan, D. J., Thornber, C., Walshaw, R., Lynn, K. J. & Trusdell, F. A. (2023a). Dynamics of magma mixing and magma mobilisation beneath Mauna Loa—insights from the 1950 AD Southwest Rift Zone eruption. *Bulletin of Volcanology* **85**(12), 1–21. https://doi.org/10.1007/S00445-023-01680-X.

Kahl, M., Mutch, E. J. F., Maclennan, J., Morgan, D. J., Couperthwaite, F., Bali, E., Thordarson, T., Guðfinnsson, G. H., Walshaw, R., Buisman, I., Buhre, S., van der Meer, Q. H. A., Caracciolo, A., Marshall, E. W., Rasmussen, M. B., Gallagher, C. R., Moreland, W. M., Höskuldsson, Á. & Askew, R. A. (2023b). Deep magma mobilization years before the 2021 CE Fagradalsfjall eruption, Iceland. *Geology* **51**(2), 184–188. https://doi.org/10.1130/G50340.1.

Kelley, K. A. & Cottrell, E. (2009). Water and the oxidation state of subduction zone magmas. *Science* **325**(5940), 605–607. https://doi.org/10.1126/science.1174156.

Kelman, M. C. (2005). *Glaciovolcanism at the Mount Cayley Volcanic Field, Garibaldi Volcanic Belt, Southwestern British Columbia*. [PhD Thesis. University of British Columbia, Vancouver (BC, Canada)]. http://hdl.handle.net/2429/16895

Kelman, M. C. & Wilson, A. M. (2024). Assessing the relative threats from Canadian volcanoes. *Canadian Journal of Earth Sciences* **61**(3), 408–430. https://doi.org/10.1139/cjes-2023-0074.

Kent, A. J., Blundy, J., Cashman, K. V., Cooper, K. M., Donnelly, C., Pallister, J. S., Reagan, M., Rowe, M. C. & Thornber, C. R. (2007). Vapor transfer prior to the October 2004 eruption of Mount St. Helens, Washington. *Geology* **35**(3), 231–234. https://doi.org/10.1130/G22809A.1.

Kent, A. J., Darr, C., Koleszar, A. M., Salisbury, M. J. & Cooper, K. M. (2010). Preferential eruption of andesitic magmas through recharge filtering. *Nature Geoscience* **3**(9), 631–636. https://doi.org/10.1038/ngeo924.

Koleszar, A. M., Kent, A. J., Wallace, P. J. & Scott, W. E. (2012). Controls on long-term low explosivity at andesitic arc volcanoes: insights from Mount Hood, Oregon. *Journal of Volcanology and Geothermal Research* **219-220**, 1–14. https://doi.org/10.1016/j.jvolgeores.2012.01.003.

Le Bas, M. J., Le Maitre, R. W., Streckeisen, A. & Zanettin, B. (1986). A chemical classification of volcanic rocks based on the total alkali-silica diagram. *Journal of Petrology* **27**(3), 745–750. https://doi.org/10.1093/petrology/27.3.745.

Leeman, W. P. & Smith, D. R. (2018). The role of magma mixing, identification of mafic magma inputs, and structure of the underlying magmatic system at Mount St. Helens. *American Mineralogist* **103**(12), 1925–1944. https://doi.org/10.2138/am-2018-6555.

Lloyd, G. E., Ferguson, C. C. & Law, R. D. (1987). Discriminatory petrofabric analysis of quartz rocks using SEM electron channelling. *Tectonophysics* **135**(1–3), 243–249. https://doi.org/10.1016/0040-1951(87)90165-X.

Mangler, M. F., Petrone, C. M., Hill, S., Delgado-Granados, H. & Prytulak, J. (2020). A pyroxenic view on magma hybridization and crystallization at Popocatépetl Volcano, Mexico. *Frontiers in Earth Science* **8**, 552516. https://doi.org/10.3389/feart.2020.00362.

Marsh, B. D. (1996). Solidification fronts and magmatic evolution. *Mineralogical Magazine* **60**(398), 5–40. https://doi.org/10.1180/minmag.1996.060.398.03.

Martin, V. M., Morgan, D. J., Jerram, D. A., Caddick, M. J., Prior, D. J. & Davidson, J. P. (2008). Bang! Month-scale eruption triggering at Santorini volcano. *Science* **321**(5893), 1178. https://doi.org/10.1126/science.1159584.

Mathews, W. H. (1958a). Geology of the Mount Garibaldi map-area, southwestern British Columbia, Canada. Part I: igneous and metamorphic rocks. *Geological Society of America Bulletin* **69**(2), 162–178. https://doi.org/10.1130/0016-7606(1958)69[161:GOTMGM]2.0.CO;2.

Mathews, W. H. (1958b). Geology of the Mount Garibaldi map-area, Southwestern British Columbia, Canada. Part II: geomorphology and quaternary volcanic rocks. *Geological Society of America Bulletin* **69**(2), 179–198. https://doi.org/10.1130/0016-7606(1958)69[179:GOTMGM]2.0.CO;2.

Menounos, B., Osborn, G., Clague, J. J. & Luckman, B. H. (2009). Latest Pleistocene and Holocene glacier fluctuations in western Canada. *Quaternary Science Reviews* **28**(21–22), 2049–2074. https://doi.org/10.1016/j.quascirev.2008.10.018.

Merrill, R. J., Bostock, M. G., Peacock, S. M., Schaeffer, A. J., & Roecker, S. W. (2022). Complex structure in the Nootka Fault Zone revealed by double-difference tomography and a new earthquake catalog. *Geochemistry, Geophysics, Geosystems*, **23**(2), e2021GC010205. https://doi.org/https://doi.org/10.1029/2021GC010205

Metcalfe, A., Moune, S., Komorowski, J. C., Kilgour, G., Jessop, D. E., Moretti, R. & Legendre, Y. (2021). Magmatic processes at La
Soufrière de Guadeloupe: insights from crystal studies and dif-
fusion timescales for eruption onset. *Frontiers in Earth Science* **9**,
617294. https://doi.org/10.3389/feart.2021.617294.

Michol, K. A., Russell, J. K. & Andrews, G. D. M. (2008). Welded block and ash flow deposits from Mount Meager, British Columbia, Canada. *Journal of Volcanology and Geothermal Research* **169**(3), 121–144. https://doi.org/10.1016/j.jvolgeores.2007.08.010.

Milman-Barris, M. S., Beckett, J. R., Baker, M. B., Hofmann, A. E., Morgan, Z., Crowley, M. R., Vielzeuf, D. & Stolper, E. (2008). Zoning of phosphorus in igneous olivine. *Contributions to Mineralogy and Petrology* **155**(6), 739–765. https://doi.org/10.1007/S00410-007-0268-7.

Morgan, D. J., Blake, S., Rogers, N. W., DeVivo, B., Rolandi, G., Macdonald, R. & Hawkesworth, C. J. (2004). Time scales of crystal residence and magma chamber volume from modelling of diffusion profiles in phenocrysts: Vesuvius 1944. *Earth and Planetary Science Letters* **222**(3–4), 933–946. https://doi.org/10.1016/j.epsl.2004.03.030.

Mourey, A. J., Shea, T. & Hammer, J. E. (2022). Preservation of magma recharge signatures in Kîlauea olivine during protracted storage. *Journal of Geophysical Research: Solid Earth* **128**(1). https://doi.org/10.1029/2022JB025523.

Mourey, A. J., Shea, T., Costa, F., Shiro, B. & Longman, R. J. (2023). Years of magma intrusion primed Kîlauea Volcano (Hawai'i) for the 2018 eruption: evidence from olivine diffusion chronometry and monitoring data. *Bulletin of Volcanology* **85**(3), 1–18. https://doi.org/10.1007/S00445-023-01633-4.

Mullen, E. K. & Weis, D. (2013). Sr-Nd-Hf-Pb isotope and trace element evidence for the origin of alkalic basalts in the Garibaldi Belt, northern Cascade arc. *Geochemistry, Geophysics, Geosystems* **14**(8), 3126–3155. https://doi.org/10.1002/ggge.20191.

Mullen, E. K. & Weis, D. (2015). Evidence for trench-parallel mantle flow in the northern Cascade arc from basalt geochemistry. *Earth and Planetary Science Letters* **414**, 100–107. https://doi.org/10.1016/j.epsl.2015.01.010.

Mullen, E. K., Weis, D., Marsh, N. B. & Martindale, M. (2017). Primitive arc magma diversity: new geochemical insights in the Cascade Arc. *Chemical Geology* **448**, 43–70. https://doi.org/10.1016/j.chemgeo.2016.11.006.

Noiriel, C., Oursin, M. & Daval, D. (2020). Examination of crystal dissolution in 3D: a way to reconcile dissolution rates in the laboratory? *Geochimica et Cosmochimica Acta* **273**, 1–25. https://doi.org/10.1016/j.GCA.2020.01.003.

Oxford Instruments (2023). *Electron Backscatter Diffraction (EBSD) Geometry*. https://www.ebsd.com/hints-and-tips/ebsd-geometry (accessed 21 November 2023)

Pankhurst, M. J., Morgan, D. J., Thordarson, T. & Loughlin, S. C. (2018). Magmatic crystal records in time, space, and process, causatively linked with volcanic unrest. *Earth and Planetary Science Letters* **493**, 231–241. https://doi.org/10.1016/j.epsl.2018.04.025.

Petrone, C. M., Bugatti, G., Braschi, E. & Tommasini, S. (2016). Pre-eruptive magmatic processes re-timed using a non-isothermal approach to magma chamber dynamics. *Nature Communications* **7**(1), 1–11. https://doi.org/10.1038/NGOMMS12946.

Phillips, M. & Till, C. B. (2022). Crustal storage and ascent history of the Mt. Shasta primitive magnesian andesite: implications for arc magma crustal flux rates. *Contributions to Mineralogy and Petrology* **177**(1), 1–27. https://doi.org/10.1007/S00410-021-01853-X.

Prior, D. J., Trimby, P. W., Weber, U. D. & Dingley, D. J. (1996). Orientation contrast imaging of microstructures in rocks using fores scatter detectors in the scanning electron microscope. *Mineralogical Magazine* **60**(403), 859–869. https://doi.org/10.1180/minmag.1996.060.403.01.

Prior, D. J., Boyle, A. P., Brenker, F., Cheadle, M. C., Austin, D., Lopez, G., Peruzzo, L., Potts, G. J., Reddy, S., Spiess, R., Timms, N. E., Trimby, P., Wheeler, J. & Zetterström, L. (1999). The application of electron backscatter diffraction and orientation contrast imaging in the SEM to textural problems in rocks. *American Mineralogist* **84**(11-12), 1741–1759. https://doi.org/10.2138/AM-1999-11-204.

Putirka, K. D. (2008). Thermometers and barometers for volcanic systems. *Reviews in Mineralogy and Geochemistry* **69**, 61–120. https://doi.org/10.2138/rmg.2008.69.3.

Putirka, K. D., Perfit, M., Ryerson, F. J. & Jackson, M. G. (2007). Ambient and excess mantle temperatures, olivine thermometry, and active vs. passive upwelling. *Chemical Geology* **241**(3–4), 177–206. https://doi.org/10.1016/j.chemgeo.2007.01.014.

Rae, A. S. P., Edmonds, M., MacLennan, J., Morgan, D., Houghton, B., Hartley, M. E. & Sides, I. (2016). Time scales of magma transport and mixing at Kîlauea Volcano, Hawai'i. *Geology* **44**(6), 463–466. https://doi.org/10.1130/G37800.1.

Rasmussen, D. J., Plank, T. A., Roman, D. C., Power, J. A., Bodnar, R. J. & Hauri, E. H. (2018). When does eruption run-up begin? Multidisciplinary insight from the 1999 eruption of Shishaldin volcano. *Earth and Planetary Science Letters* **486**, 1–14. https://doi.org/10.1016/j.epsl.2018.01.001.

Read, P. B. (1977a). Geology of Meager Creek geothermal area, British Columbia. *Geological Survey of Canada Open File* **603**.

Read, P. B. (1977b). Meager Creek volcanic complex, southwestern British Columbia. *Geological Survey of Canada Paper* **77-1A**, 277–281.

Read, P. B. (1990). Mount Meager Complex, Garibaldi Belt, South-Western British Columbia. *Geoscience Canada* **17**(3), 167–174. https://journals.lib.unb.ca/index.php/GC/article/view/3672.

Reed, S. J. B. (2005) *Electron microprobe analysis and scanning electron microscopy in geology*. Cambridge (United Kingdom): Cambridge University Press.

Reyes, A. V. & Clague, J. J. (2004). Stratigraphic evidence for multiple Holocene advances of Lillooet Glacier, southern Coast Mountains, British Columbia. *Canadian Journal of Earth Sciences* **41**(8), 903–918. https://doi.org/10.1139/e04-039.

Rhodes, J. M., Dungan, M. A., Blanchard, D. P. & Long, P. E. (1979). Magma mixing at mid-ocean ridges: evidence from basalts drilled near 22° N on the Mid-Atlantic Ridge. *Tectonophysics* **55**(1-2), 35–61. https://doi.org/10.1016/0040-1951(79)90334-2.

Riddihough, R. (1984). Recent movements of the Juan de Fuca plate system. *Journal of Geophysical Research* **89**(B8), 6980–6994. https://doi.org/10.1029/JB089iB08p06980.

Roberti, G., Ward, B., van Wyk de Vries, B., Friele, P., Perotti, L., Clague, J. J. & Giardino, M. (2018). Precursory slope distress prior to the 2010 Mount Meager landslide, British Columbia. *Landslides* **15**(4), 637–647. https://doi.org/10.1007/s10346-017-0901-0.

Roeder, P. L. & Emslie, R. F. (1970). Olivine-liquid equilibrium. *Contributions to Mineralogy and Petrology* **29**(4), 275–289. https://doi.org/10.1007/BF00371276.

Russell, J. K., Stewart, M., Wilson, A. & Williams-Jones, G. (2021). Erup-
tion of Mount Meager, British Columbia, during the early Fraser
glaciation. *Canadian Journal of Earth Sciences* **58**(10), 1146–1154.
https://doi.org/10.1139/cjes-2021-0023.

Russell, J. K., Edwards, B. R., Williams-Jones, G. & Hickson, C. (2023). Pleistocene to Holocene volcanism in the Canadian Cordillera. *Canadian Journal of Earth Sciences* **60**, 1443–1466. https://doi.org/10.1139/CJES-2023-0065.

Ruth, D. C. S., Costa, F., Bouvet De Maisonneuve, C., Franco, L., Cortés, J. A. & Calder, E. S. (2018). Crystal and melt inclusion timescales reveal the evolution of magma migration before erup-
tion. *Nature Communications* **9**(1), 2657. https://doi.org/10.1038/
s41467-018-05086-8.
Saunders, K., Blundy, J., Dohmen, R. & Cashman, K. (2012). Linking petrology and seismology at an active volcano. *Science* **336**(6084), 1023–1027. https://doi.org/10.1126/science.1220066.

Schneider, C. A., Rasband, W. S. & Eliceiri, K. W. (2012). NIH image to ImageJ: 25 years of image analysis. *Nature Methods* **9**(7), 671–675. https://doi.org/10.1038/nmeth.2089.

Shea, T., Costa, F., Krimer, D. & Hammer, J. E. (2015). Accuracy of timescales retrieved from diffusion modeling in olivine: a 3D perspective. *American Mineralogist* **100**(10), 2026–2042. https://doi.org/10.2138/am-2015-5163.

Stasiuk, M. V., Russell, J. K. & Hickson, C. J. (1996). Distribution, nature, and origins of the 2400 BP eruption products of Mount Meager, British Columbia: linkages between magma chemistry and eruption behaviour. *Geological Survey of Canada Bulletin* **486**, 27. https://doi.org/10.4095/208165.

Stockstill, K. R., Vogel, T. A. & Sisson, T. W. (2003). Origin and emplacement of the andesite of Burroughs Mountain, a zoned, large-volume lava flow at Mount Rainier, Washington, USA. *Journal of Volcanology and Geothermal Research* **119**(1–4), 275–296. https://doi.org/10.1016/S0377-0273(02)00358-X.

Thorkelson, D. J., Madsen, J. K. & Sluggett, C. L. (2011). Mantle flow through the northern cordilleran slab window revealed by volcanic geochemistry. *Geology* **39**(3), 267–270. https://doi.org/10.1130/G31522.1.

Thornber, C. R. & Huebner, S. J. (1985). Dissolution of olivine in basaltic liquids: experimental observations and applications. *American Mineralogist* **70**(9–10), 934–945.

Toplis, M. J. (2005). The thermodynamics of iron and magnesium partitioning between olivine and liquid: criteria for assessing and predicting equilibrium in natural and experimental systems. *Contributions to Mineralogy and Petrology* **149**, 22–39. https://doi.org/10.1007/s00410-004-0629-4.

Ubide, T., Mollo, S., Zhao, J. X., Nazzari, M. & Scarlato, P. (2019). Sector-zoned clinopyroxene as a recorder of magma history, eruption triggers, and ascent rates. *Geochimica et Cosmochimica Acta* **251**, 265–283. https://doi.org/10.1016/j.gca.2019.02.021.

Unnsteinsson, T., Flowers, G. & Williams-Jones, G. (2024). Formation and persistence of glaciovolcanic voids explored with analytical and numerical models. *Journal of Glaciology* **1-15**, 1–15. https://doi.org/10.1017/jog.2024.8.

Venugopal, S., Moune, S., Williams-Jones, G., Druitt, T., Vigouroux, N., Wilson, A. & Russell, J. K. (2020). Two distinct mantle

sources beneath the Garibaldi Volcanic Belt: insight from olivine-hosted melt inclusions. *Chemical Geology* **532**, 119346. https://doi.org/10.1016/j.chemgeo.2019.119346.

Voloschina, M., Métrich, N., Bertagnini, A., Marianelli, P., Aiuppa, A., Ripepe, M. & Pistolesi, M. (2023). Explosive eruptions at Stromboli volcano (Italy): a comprehensive geochemical view on magma sources and intensity range. *Bulletin of Volcanology* **85**(6), 34. https://doi.org/10.1007/s00445-023-01647-y.

Walker, L. (2003). *Late Quaternary Glacier Fluctuations at Lillooet, Diadem, and Berendon Glaciers, Coast Mountains, British Columbia* [PhD Thesis. Simon Fraser University, Canada]. https://summit.sfu.ca/item/8591

Walowski, K. J., Wallace, P. J., Cashman, K. V., Marks, J. K., Clynne, M. A. & Ruprecht, P. (2019). Understanding melt evolution and eruption dynamics of the 1666 C.E. Eruption of cinder cone, Lassen Volcanic National Park, California: insights from olivine-hosted melt inclusions. *Journal of Volcanology and Geothermal Research* **387**, 106665. https://doi.org/10.1016/j.jvolgeores.2019.106665.

Welsch, B., Faure, F., Famin, V., Baronnet, A. & Bachèlery, P. (2013). Dendritic crystallization: a single process for all the textures of olivine in basalts? *Journal of Petrology* **54**(3), 539–574. https://doi.org/10.1093/petrology/egs077.

Welsch, B., Hammer, J. & Hellebrand, E. (2014). Phosphorus zoning reveals dendritic architecture of olivine. *Geology* **42**(10), 867–870. https://doi.org/10.1130/G35691.1.

Wieser, P., Petrelli, M., Lubbers, J., Wieser, E., Ozaydin, S., Kent, A. & Till, C. (2022). Thermobar: an open-source Python3 tool for thermobarometry and hygrometry. *Volcanica* **5**(2), 349–384. https://doi.org/10.30909/vol.05.02.349384.

Wilson, A. M. & Russell, J. K. (2017). Lillooet glacier basalts, southwestern British Columbia, Canada: products of quaternary glaciovolcanism. *Canadian Journal of Earth Sciences* **54**(6), 639–653. https://doi.org/10.1139/cjes-2016-0201.

Wilson, A. M. & Russell, J. K. (2018) Quaternary glaciovolcanism in the Canadian Cascade volcanic arc—Paleoenvironmental implications. In: Poland M., Garcia M., Camp V. & Grunder A. (eds) *Field Volcanology: A Tribute to the Distinguished Career of Don Swanson*, Special Paper 538. Boulder, Colorado, US:Geological Society of America, pp.133–157. https://doi.org/10.1130/2018.2538(06).

Zellmer, G. F. (2021). Gaining acuity on crystal terminology in volcanic rocks. *Bulletin of Volcanology* **83**(11), 1–8. https://doi.org/10.1007/s00445-021-01505-9.
<figure><img src="image_16.png" /></figure>

Deltech Furnaces

Sustained operating
temperatures to 1800°
Celsius

www.deltechfurnaces.com

<figure><img src="image_17.png" /><figcaption>Gas Mixing System</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>An ISO 9001:2015 certified company</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>ASME NQA-1 2008 Nuclear Quality Assurance</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Control systems are certified by Intertek UL508A compliant</figcaption></figure>

<caption>Bottom Loading Vertical Tube</caption>