Cheakamus basalt lavas, British Columbia: a Pleistocene record of rapid, continuous eruption within a mountainous drainage system

The Cascade volcanic arc comprises a chain of volcanoes spanning ∼1250 km from northern California, USA, to southwest British Columbia, Canada. The Garibaldi Volcanic Belt (GVB) is the northern extent of the Cascade volcanic arc (Fig. 1, inset). Whereas the southern Cascades feature large stratovolcanoes in Washington, Oregon, and California (i.e., Mt. St Helens and Mt. Rainier) reflecting more focused volcanism (Hickson 1994; Hildreth 2007; Ramsey and Siebert 2017), the GVB is characterized by complex and discreet volcanic fields (i.e., Mt. Meager Volcanic Complex and Mt. Cayley Volcanic Field). The GVB differs further in that it overlies crystalline bedrock of the Coast Mountain Range where local relief from valley floors to summits is commonly 1500–2000 m. A history of repeated glaciations has left a steep, mountainous, and highly dissected landscape of fjords, glacial cirques, and hanging- and U-shaped valleys.

Large-volume, high-standing, partially dissected stratovolcanoes dominate the GVB; however, smaller monogenetic basaltic volcanoes and lavas are numerous (Wood and Kienle 1990; Hickson 1994; Hildreth 2007). The effusive eruptions of these low-viscosity basaltic lavas represent geological events with the capacity to transform landscapes over very short timescales. These eruptions and their deposits can result in rapid and major derangement of landscapes and can control the subsequent recovery and evolution of that landscape. This is especially true in the GVB where the profound topographic relief can channel lava into major valley drainages, causing blockages and derangement of active and rapidly evolving glacial, peri-glacial, and fluvial systems.

The Cheakamus lavas, first described by Mathews (1948, 1958), are one such example of an effusive volcanic eruption exploiting and modifying drainage systems in a mountainous terrain, as well as impacting the post-glacial recovery of the landscape. Here, we present a geological map for the distribution of Cheakamus lavas that includes an inferred vent location, an estimate of the original surface area inundated by lava, and the total volume erupted. Our map distinguishes three separate phases of eruption based on stratigraphic relationships and minor variations in petrographic and geochemical properties. Our field work includes sampling for paleomagnetic study and ⁴⁰Ar/³⁹Ar dating with the aim of constraining both the duration and age of eruption. We suggest that the Cheakamus lavas were erupted continuously over a short period of time during the early stages of the Fraser glaciation (∼26–18 ka; Clague and Ward 2011) when the major drainages were largely ice-free.

The Cheakamus lavas are situated immediately east of the Mt. Cayley Volcanic Field (Kelman et al. 2002) and north of Mt. Garibaldi (Green et al. 1988), both of which are dominated by calc-alkaline intermediate to felsic volcanic rocks. Volcanic edifices in this part of the GVB span an eruptive history from 2.7 Ma (Mt. Price; Mathews 1948, 1958; Green et al. 1988) to 10 ka (Ring Creek lava; Brooks and Friele 1992). Slag Hill and Ring Mountain are intermediate calc-alkaline edifices within the Mt. Cayley Volcanic Field and are the closest volcanoes to the inferred vent of the Cheakamus basalts (Fig. 1).

The Cheakamus lavas are some of the youngest Quaternary lavas within the GVB (Mathews 1948, 1958; Green 1990; Wilson and Russell 2018) and mainly occupy the present-day Callaghan and Cheakamus River valleys between the towns of Squamish and Whistler, British Columbia. The lavas stretch ∼31 km from their source near Conflict Lake at the top of the Callaghan Valley to their present-day terminus near Daisy Lake (Fig. 1). They have individual flow thicknesses of ∼2–15 m and several exposures comprise multiple (N = 5–7) lavas having combined thicknesses of ∼80 m (Fig. 2). Extensive Quaternary and Holocene glacial sediment and colluvium cover the tributary valleys upstream of the lavas occupying the Cheakamus River valley. The lava surfaces within the Cheakamus River valley are largely till-free. Extensive till cover within the Callaghan Valley obscures most of the Cheakamus lavas (Blaise-Stevens 2008); exceptions include exposures in the banks of Callaghan Creek, the shorelines of Callaghan and Conflict lakes, and a ∼45 m section at Alexander Falls comprising multiple individual lavas (Fig. 2).

Despite their proximity to a major transportation corridor and the municipality of Whistler, the age, distribution, and eruptive history of the Cheakamus lavas are not well constrained. The lavas were first described by Mathews (1948), who named them, determined their lithology (i.e., olivine (Ol) basalt), mapped their distribution, and postulated a source at Callaghan Lake. He estimated an eruptive volume of 1.25 km³. Mathews (1948, 1958) used outcrop morphology, the presence of interbedded sediments, and the extent of dissection to suggest a Late Pleistocene to Recent age and postulated a protracted eruptive history spanning waxing and waning of the Cordilleran ice sheet (i.e., Fraser glaciation). He also proposed a glaciovolcanic origin for the stratigraphically highest lavas that he noted as having an “esker-like” morphology, abundant enigmatic radially oriented columnar jointing, and partially palagonitized basal contacts in some exposures. Mathews (1958) interpreted these high-standing outcrops of basalt as lavas that exploited subglacial drainage systems at the close of the Fraser glaciation.

Subsequent work by Green (1977) provided more detailed petrological, geochemical, and isotopic characterization of the Cheakamus lavas (Green et al. 1988; Green 1990). Green (1977, 1981) divided Mathews’ Cheakamus River valley basalts into multiple stages based on some of the same criteria as Mathews (1948, 1958) combined with minor variations in geochemical composition and phenocryst abundance. He proposed four main phases of eruption and ascribed chemical variations to fractionation processes involving Ol, clinopyroxene (Cpx), and plagioclase (Pl) (Green 1977). Green (1977, 1981, 1988) also reported on two radiometric age estimates for the Cheakamus lavas. A single sample of lava collected near Daisy Lake yielded an imprecise K–Ar age of 50 ± 50 ka. A ¹⁴C age estimate based on an organic sample collected from sediment situated between two lavas returned an age of 34 000 ± 800 years (Green 1981; McNeely 1989). However, the authors noted the low amount of ¹⁴C and the uncertain origin of the fine organic detritus collected (McNeely 1989). Indeed, if the material that was dated was not formed in situ but washed in from one or more distal sources (e.g., Punning and Rajamäe 1993), the ¹⁴C date would represent a maximum age for the upper flow, leaving the age of the lower lava as unconstrained. Our attempts to redate the same sedimentary layer by ¹⁴C returned only modern ages. Green (1977, 1988) accepted the ¹⁴C age estimate and ascribed the interbedded sediments to the Olympia interstade (∼50–25 ka; Clague and Ward 2011) and the underlying lavas to pre-Salmon Springs glaciation (considered ∼70 ka at the time of the study; Easterbrook et al. 1981).

The Cheakamus lavas occupy the present-day Conflict Lake drainage and Callaghan and Cheakamus River valleys (Fig. 1). There is no discernible vent, but the distribution of lava suggests a source located near Conflict Lake at an elevation of ∼1380 masl. The present-day drainage is occupied by a highly glaciated plateau of lava extending from Conflict Lake into the Callaghan River valley. Rare exposures in the Callaghan drainage suggest the paleo-valley is filled by stacks of lavas reaching thicknesses of ≥45 m (Fig. 2C; Table 1) and extending ∼14 km downstream before reaching the Cheakamus Valley. From the junction of the Callaghan and Cheakamus River valleys, at an elevation of 506 masl, the lavas fill the Cheakamus River valley to a maximum thickness of ∼80 m as observed at Brandywine Falls (Fig. 2D; Table 1) and extend downstream another ∼12 km before terminating south of Daisy Lake (370 masl). The total elevation drop from the vent area to the lava terminus is ∼1.1 km.

In several instances, where two drainage systems intersect (e.g., Conflict Lake–Callaghan and Callaghan–Cheakamus), lava flowing into the new drainage has backed up in the upstream direction, while the main flowage has been downslope (Fig. 1). In the Callaghan Valley, bedrock constrictions narrow the valley and the subsequent width of the lavas therein. The Cheakamus lavas cover a total area of ∼28 km², and based on paleo-valley cross-section reconstructions and thickness measurements at well-exposed sections, the basalts have an estimated present-day volume of 1.45 km³. We estimate the original, pre-erosional volume of Cheakamus basalt lavas to be ∼1.65 km³ (Table 1).

The Cheakamus lavas are Ol and Pl porphyritic basalts exhibiting little petrographic variation other than in (micro-)phenocryst content (14%–25%). A number of exposures comprise multiple individual lavas sometimes featuring interbeds or discontinuous lenses of clastic sediments (Figs. 2A and 3). Geochemical variation between lavas is limited (see below) and, thus, our stratigraphic framework is based mainly on field relationships (i.e., interbeds of sediments), geographic distributions (i.e., drainage systems), and morphological and surface features of lavas (e.g., glacial scouring; Fig. 2B). We define four mappable units within the Cheakamus basalts (Figs. 1 and 3; Table 2), including Early Brandywine phase (EBp) lavas, interflow sediments, Brandywine Falls phase (BFp) lavas, and Upper Callaghan phase (UCp) lavas. The volcanic successions preserve a number of primary textures and features that inform on the syn-eruptive paleo-environments and processes, including high-relief, ice-free drainages, drainages occupied by substantial river systems, and, potentially, localized masses of valley-filling glaciers.

The EBp lavas are best exposed in the cliffs at Brandywine Falls (Fig. 2D) and in escarpments surrounding Daisy Lake. The EBp lavas occupy the lowest stratigraphic position in the Cheakamus and Callaghan valleys and rest on bedrock in multiple places—notably at the base of Brandywine Falls and on Daisy Lake Island. The EBp comprises 2–6 individual basalt lavas that form a broad subaerial flow field filling the paleo-Cheakamus River valley to a maximum observed thickness of ∼80 m (Fig. 2D). The majority of lavas have smooth pahoehoe surfaces and feature 0.5–1 m of basal autobreccias. The EBp lavas have well-developed lower colonnades below entablatures of equal thickness (Figs. 2B and 2D) and, where exposed in the modern-day Brandywine and Cheakamus River valleys, are seen to have poorly formed or entirely lack an upper colonnade (Fig. 2). Thick sections of EBp lavas are highly eroded by later events, especially in the southernmost exposures around Daisy Lake. The effect of this late erosion is to form high-standing ridges and plateaus that are incised by a series of channels. In some cases, dissected surfaces of EBp lavas are heavily scoured but always lack glacial sediment cover (Fig. 3; Table 2).

Discontinuous beds and lenses of clastic sediments are found interlayered between basalt lavas in outcrops exposed in the modern-day Cheakamus River valley. In the northern and southern reaches of the Cheakamus Valley (Figs. 1 and 3), the interflow sediment is continuously visible over scales of ∼10 m and reaches thicknesses of 75 cm (Fig. 2A). In the Brandywine Falls area, sediment presents as discontinuous 1–2 m lenses ∼40 cm thick. Interflow sediments are notably absent from lava sequences exposed in the Callaghan Valley.

The beds of interflow sediments vary in character. In places, the unit comprises weakly stratified medium-fine sand grading normally to fine sand and silt, with discontinuous lenses of well-rounded pebble to cobble gravel. Some cobbles have a well-developed flat-iron shape, with striated facets that suggest these cobbles were previously glacially transported (i.e., Martini et al. 2011). Variably glassy and palagonitized bases of overlying lavas (Fig. 2A; Wilson and Russell 2018) suggest that the sediment was wet when overrun by subsequent lavas (i.e., Fuller 1931). In the southern reaches of the Cheakamus Valley, the interflow sediment is characterized by discreet 1–2 m pockets of silt up to 40 cm thick, often mixed with clasts of basal breccia showing glassy rinds with minor palagonitization and sediment baking (Fig. 2A). Where interflow sediments are absent, the contact between lavas commonly features a 0.2–1.5 m thick basal autobreccia (i.e., related to the upper lava).

The interflow sediments suggest some amount of time, no matter how short, between emplacement of the lower (i.e., earlier) and upper (i.e., later) lavas. We use the interflow sediments to separate the lavas representing the early phase of eruption (EBp) from later lavas that continue to fill the paleo-Cheakamus River valley that we define as the BFp (see below; Table 2).

BFp lavas occupy the highest stratigraphic position in the Cheakamus and lower Callaghan valleys, and conformably overlie discontinuous layers and lenses of interflow sediment. Locally, basal contacts of the BFp lavas commonly feature bulbous, glassy 5–10 cm protrusions into the underlying sediment and can be partially palagonitized (Fig. 2A). As described above, BFp lavas typically have 0.2–1.5 m thick basal autobreccias where underlying interflow sediment is absent. In the southern Cheakamus Valley, the BFp comprises a single ∼3–8 m thick, subaerial basalt lava. That lava has similar features as found in EBp lavas, including a lower colonnade, a well-developed entablature, and no upper colonnade (Fig. 2B). It also has an eroded upper surface and, unlike lavas in the Callaghan Valley, is not blanketed by till (Blaise-Stevens 2008). To the north, BFp lavas enter the Cheakamus Valley from the tributary Callaghan Valley and travel ∼1.2 km upstream to an elevation of 543 m (30 m above the junction of the two valleys), marking the Northern Terminus of Cheakamus lavas (Figs. 1 and 3B). There, the exposure comprises two interfingering, 2–5 m thick flows. At this location, lavas display fine, 4–10 cm cube jointing and 20–50 cm curving joints and form tall, narrow, discontinuous anastomosing ridges with scoured surfaces (cf. Mathews 1948, 1958). Lavas are perched ∼35 m above the valley floor on the western valley wall but are absent on the eastern valley wall at equivalent elevations and absent in the present-day valley bottom.

The northernmost reaches of the Callaghan Valley, as well as the Conflict Lake drainage, are filled by 0.5–2 m thick deposits of till and ice-contact deposits of sand and gravel (Blaise-Stevens 2008); exposures of the underlying volcanic rocks are thus limited. Rare exposures suggest that the majority of the Callaghan Valley is filled by BFp and EBp lavas (Fig. 2C); however, there is no evidence of interflow sediment between the lavas (as seen in lavas within the Cheakamus Valley; Table 2). The EBp lavas are themselves overlain by lavas of the UCp that appear to source from Conflict Lake and persist for ∼5 km before being lost to till cover in the Callaghan Valley. Where exposed, they form an ∼4–8 m thick, light grey basaltic andesite lava featuring coarse 1–2 m columns and a heavily glaciated surface (Table 2). The contact between UCp lavas and the underlying EBp/BFp lavas is not exposed. The UCp is the most proximal phase to the vent location and most likely the last eruptive phase; however, poor exposure and intervening distances obscure the UCp’s stratigraphic relationship to lavas in the Cheakamus Valley (i.e., EBp).

The Cheakamus lavas are Ol and Pl porphyritic, and show only minor petrographic variations. Most samples are characterized by variable abundances of Ol (∼0.25–1 mm) and Pl (∼0.5–2 mm) microphenocrysts (or glomerocrysts) and tabular, euhedral to subhedral, sieve-textured Pl xenocrysts (or possibly antecrysts; 2–6 mm; 3%–7%; Table 3; Fig. 4) within a holocrystalline groundmass of Pl, Cpx, and Fe–Ti oxides. Some samples contain groundmass Ol (Figs. 4A and 4B), while others lack groundmass Ol altogether and are reported as containing pigeonite (Nicholls et al. 1982). Several lavas exposed in the lower parts of the Callaghan River valley feature ophitic masses of pleochroic titanaugite within a weakly trachytic groundmass (Fig. 4D). The EBp and BFp lavas are most similar petrographically, although groundmass Ol is slightly more common in EBp lava (Fig. 4A; Table 3). The UCp lavas show the greatest difference in that Pl phenocrysts are more abundant (∼6% compared to ∼2% in EBp lavas; Fig. 4C). UCp lavas also have higher abundances of xenocrystic (or antecrystic) Pl and glomerophyric clots of Pl and Ol are more common than in other phases (Fig. 4C; Table 3). The groundmass is often trachytic, and groundmass Ol is rare to absent.

Whole rock major and trace element geochemical compositions were measured for 32 samples by X-ray fluorescence analysis of fused discs and by inductively coupled plasma-mass spectroscopy, respectively (see Appendix A). The sample suite is representative of all phases of the Cheakamus lavas (Fig. 1; Table 2).

All samples plot as subalkaline basalt (EBp and BFp phases) to basaltic andesite (UCp) (Fig. 5A; Irvine and Baragar 1971; Le Bas et al. 1986). Chemically, all of the lavas are hypersthene normative and have 8%–22% normative Ol. The main chemical variations are expressed by SiO₂ (∼49%–53 wt.%), MgO (6–8 wt.%), and FeO (9.5–11.8 wt.%) and there is strong linear negative covariation between MgO and FeO with SiO₂. The chemical data show an ordered progression to lower values of MgO + FeO from EBp lavas to BFp lavas and to the most SiO₂-rich UCp lavas (Fig. 5B). The geochemical trend in CaO contents is more complicated. EBp lavas have the highest CaO and MgO + FeO contents and show a linear decrease in both that is continuous with the BFp lavas; the UCp lavas overlap some of the BCp lavas but then show a pronounced increase in CaO with decreasing MgO + FeO (Fig. 5C).

Rare earth element (REE) compositions of Cheakamus lavas are shown in Fig. 6A normalized to primitive mantle composition (Sun and McDonough 1989). The REE patterns are parallel for samples from all three phases (i.e., EBp, BFp, and UCp) and show little relative enrichment or depletion between samples suggesting a common mantle source and limited differentiation. Ratios of the incompatible elements Y and Zr are constant to within analytical error for all three phases of lavas and across the full range of MgO + FeO contents (Fig. 6B) lending further support to a common magmatic origin. Conversely, trace element ratios of Ni/Zr and Sr/Zr show pronounced variations significantly greater than analytical error that are strongly (negatively and positively, respectively) correlated to MgO + FeO (Figs. 6C and 6D). If we accept MgO + FeO as a metric of differentiation, then Figs. 5C and 5D suggest that the magmatic differentiation processes are accompanied by a depletion or loss of Ni and an enrichment or gain in Sr.

Pearce element ratios (PERs) offer a robust and effective means of testing petrologic hypotheses (Pearce 1968; Russell and Nicholls 1988; Stanley and Russell 1989; Russell et al. 1990; Nicholls and Russell 2016). Because they preserve geochemical information in absolute, stoichiometrically defined units, PERs remove the distorting effects of closure (i.e., constant-sum problem; Chayes 1960). Thus, PER diagrams can show absolute changes in geochemical composition attending magmatic processes (Russell and Nicholls 1988). Element ratio diagrams can be tailored to test for certain mineral phases that are involved in magmatic sorting (i.e., crystal fractionation or accumulation). For a PER to be effective, the denominator constituent must remain conserved by processes affecting the chemical evolution of the system (Pearce 1968).

The REE data (Fig. 6A) and ratios of incompatible trace elements (Fig. 6B) indicate that the Cheakamus lavas originated from a single magmatic system or from chemically indistinguishable magmas. In basaltic magmas, there are usually several elements (minor and trace) that are largely incompatible and likely to be conserved during most differentiation processes. For our purposes, we have chosen Ti as a conserved denominator. Ti is well measured and conserved relative to the observed phenocryst (Ol and Pl) and xenocryst (Pl) assemblages (Russell and Nicholls 1988).

Figure 7 is a PER diagram designed to test whole rock chemical analyses for the relative chemical effects of physical sorting processes involving Ol and/or Pl. Mineral response vectors show the effects of gaining or losing molar equivalents of Ol and Pl. Variation due to crystal-sorting of Ol alone will cause the data to spread vertically upwards (accumulation) or downwards (fractionation). The effects of Pl sorting are to create a horizontal trend, and combinations of Ol and Pl will create chemical trends where the slope is a quantitative representation of the relative proportions of the two phases. The additional attribute of the diagram is that the chemical effects of Cpx are null. The data are plotted with 1σ error envelopes resulting from the propagation of the analytical uncertainties through the ratio calculations. The error ellipses account for, both, the analytical uncertainty and the induced correlation caused by the ratios sharing a common denominator.

The three suites of Cheakamus lavas are colour coded as in previous diagrams, and are well ordered in terms of their stratigraphic position. Normal differentiation processes would drive data from the most primitive composition towards lower values of X and Y. Here, the chemical data behave differently, where, relative to the most Mg-rich sample (AB-20-46), the data move to higher X-axis values and lower Y-axis values. The most Mg-rich samples, belonging to the EBp lavas, plot in the upper left corner of the figure; the BFp lavas define a horizontal trend to higher X-axis values, suggesting little loss of Ol and accumulation of Pl. The most SiO₂-rich samples, belonging to the UCp lavas, define a trend to higher X-axis and lower Y-axis values. Relative to the most magnesian lavas (i.e., EBp), the UCp lavas are consistent with minor Ol fractionation and a net accumulation of Pl. We cannot rule out that BFp and UCp lavas may have lost Pl (i.e., crystal fractionation), but in both cases the amount lost must be subordinate to the total amount of Pl gained (e.g., accumulation, contamination, etc.).

We performed paleomagnetic sampling and analysis of the Cheakamus basalts to constrain the duration of volcanism. The sample sites span the entire stratigraphic and geographic range of the Cheakamus basalts within the Conflict Lake, Callaghan, and Cheakamus drainages (Figs. 1 and 3; Table 4). In the 2020 field season, we sampled a total of 17 sites (no. of cores = 128) across the field area for paleomagnetic study and, in 2021, added another 17 sites (a total of 290 individual core samples), giving us a total of 34 sites (Table 4).

Every effort was made to sample only intact outcrops, thereby avoiding areas suspected of post-emplacement movement. At each site, approximately 8 standard 2.5 cm diameter paleomagnetic cores were collected. Each sample was oriented using a magnetic compass, and where possible a sun compass was used in conjunction with magnetic orientations. In many cases, we sampled multiple lavas exposed in a single outcrop (e.g., see Lucille Lake, Northern Terminus; Table 4).

All samples were analysed at the paleomagnetic laboratory at the University of Lethbridge, Alberta. Magnetic susceptibility was determined with a Sapphire Instruments (SI-2B) susceptibility meter. The magnetization of each sample was measured with an AGICO JR-6A spinner magnetometer before demagnetization and again after each level of stepwise demagnetization. Cheakamus basalt samples were held in magnetic shields following field collection and between laboratory measurements. All samples were subjected to alternating field (AF) demagnetization, performed using an ASC Scientific D-2000 demagnetizer with a three-axis manual tumbler and carried out at 10 millitesla (mT) steps (up to 100 mT and in a few cases up to 200 mT). Thermal demagnetization was carried out at 100, 200, 300, 400, 500, 550, and 580 °C, using an ASC Model TD48 dual-chamber thermal demagnetizer to confirm that AF demagnetization was sufficient to resolve the primary remanence. Directions of characteristic magnetization were determined for each sample by principal component analysis (Kirschvink 1980) using Remasoft version 3.0 (Chadima and Hrouda 2006). Mean directions of characteristic magnetization were calculated for each site and an overall mean was also calculated (Table 4; Appendix C). All samples were subjected to stepwise AF demagnetization and principal component analysis and mean directions were calculated from the AF data only.

All paleomagnetic samples are normally magnetized (Fig. 8A). The mean directions of sites are consistent, where all sites plot within the uncertainty circle (α₉₅) about the overall mean (N = 34) with a declination of 345.2° and inclination of 73.0°, α₉₅ = 1.3°, and Fisher’s precision factor (k) = 319 (Fisher 1953), (Table 4; Fig. 8A). Appendix C shows the typical, well-behaved demagnetization characteristics of Cheakamus lavas for samples from EBp, BFp, and UCp.

All paleomagnetic samples produced stable directions of magnetization following stepwise demagnetization. AF demagnetization, applied to all specimens, exhibited linear decay to the origin on orthogonal projections (Appendix C). The median destructive fields range from 20 to 80 mT and final directions are often stable up to 200 mT AF demagnetization (see sample CBB017A, Appendix C), indicative of single-domain magnetite. Thermal demagnetization showed complete demagnetization at 550–580 °C, typical of fine-grained (single domain) magnetite as the magnetic carrier (see samples CBB212B and CBB351B, Appendix C). Samples from a few flows revealed softer multi-domain magnetite upon thermal demagnetization (CBB017B and CBB131B). Note that single-domain magnetite is the optimal magnetic mineral and grain size for producing reliable paleomagnetic directions.

The common mean pole direction obtained for all volcanic units, from the stratigraphically lowest to highest lavas, implies a common eruption age. Furthermore, the data from each site showed no indication of rotation or tilt as the mean direction is within the expected paleosecular variation. All site means (N = 34) overlap within a 95% confidence interval, indicating that the Cheakamus basalts were emplaced within the same paleomagnetic moment and were not erupted over a protracted period (cf. Mathews 1958; Wilson and Russell 2018). In locales where the entire stratigraphic section could be sampled (i.e., Brandywine Falls and Alexander Falls), we found no systematic drift of inclination and declination between the first and last lavas (Table 4; Fig. 8B), indicating that the eruption duration was remarkably short (<2000 years; e.g., Hagstrum and Champion 2002). This refutes previous hypotheses of a protracted eruption history (e.g., Mathews 1948; Green et al. 1988; Wilson and Russell 2018).

Previous estimates of absolute age for the Cheakamus basalts are few, have significant uncertainties, and do not adequately constrain the eruption age (Green et al. 1988; McNeely 1989). We targeted the highest and lowest stratigraphic units for ⁴⁰Ar/³⁹Ar dating using outcrops exposed at Alexander Falls (Fig. 2C) and Brandywine Falls (Fig. 2D) in the Callaghan and Cheakamus River valleys, respectively (Table 5; Appendix D). The highest stratigraphic unit at Brandywine Falls, located on the cliff 70 m above the base of the section (Fig. 3), returned a plateau age of 13.1 ± 7.2 ka (1σ uncertainty). The lowest lava in the section at Brandywine Falls returned a plateau age of 23.9 ± 7.9 ka, which overlaps the younger age estimate. The lowest exposed flow at Alexander Falls returned a plateau age of 13.4 ± 5.9 ka (Fig. 9).

The Cheakamus basalts are transitional, comagmatic Ol–Pl porphyritic basalts originating from a single magma batch (Figs. 5–7). The lack of Cpx in the phenocryst assemblage suggests the source magma was stored at shallow depths (e.g., Wilson and Russell 2017; Harris and Russell 2022). The chemical diversity of the basalts, shown via PERs (Fig. 7), indicates differential sorting of Ol and Pl during storage and ascent. The Ol population formed solely from magmatic crystallization, but the Pl population is both phenocrystic and xenocrystic in origin (Fig. 4).

The UCp lavas are somewhat petrographically and chemically distinct relative to the earlier Cheakamus basalt lavas. They have higher SiO₂ contents (i.e., basaltic andesite), greater (micro-)phenocryst abundances (∼21% Pl, ∼6% Ol), and more xenocrystic Pl. The Pl accumulation signal seen in Fig. 7 is strongest in the UCp lavas and is likely due partly to entrainment of foreign Pl. They are the latest lavas to erupt. One explanation for their petrographic and chemical difference is that the UCp lavas may represent late tapping of magma from the sidewalls of the reservoir. Magma at the margins of the reservoir would be cooler, more crystallized, and more likely to entrain xenocrystic or antecrystic material.

Field relationships indicate that the eruption was continuous and that there was no appreciable, long-term time-breaks during eruption other than local syn-eruption fluvial or glacio-fluvial sedimentation between lavas within the Cheakamus Valley. The eruption poured 1.65 km³ of low-viscosity lava into the high-relief mountainous landscape of the Coast Mountains, inundating the Callaghan and Cheakamus drainages. Modern voluminous basaltic eruptions in Hawaii and Iceland are observed to have peak and average discharge rates of 280 and ∼50 m³ s⁻¹ (1984 Mauna Loa eruption) and 350 and ∼160 m³ s⁻¹ (2014–2015 Holuhraun eruption), respectively (Pedersen et al. 2017; Plank et al. 2021). Applying these peak and average discharge rates to the Cheakamus lavas would imply corresponding eruption durations of 65–380 or 55–120 days, respectively. These durations are like those observed in recent effusive eruptions, including the 2021 Cumbre Vieja eruption (85 days; Carracedo et al. 2022), the 2018 Kilauea eruption (98 days; Neal et al. 2019), and the 2021 Fagradalsfjall eruption (183 days; Pedersen et al. 2022). In many instances, effusive eruptions of large volumes of basalt feature minor pauses of days to weeks (e.g., Wolfe et al. 1988; Mattox et al. 1993) that would potentially extend the total time involved in the emplacement of these lavas.

There are two additional and related features in the distribution of the Cheakamus lavas that pertain to their emplacement. First, there are several instances in the Callaghan Valley where basalt has inundated tributary drainages and been forced upstream to higher elevations, including (i) at Callaghan Lake, (ii) the west side of the valley between Callaghan Lake and Alexander Falls, (iii) at Alexander Falls, and (iv) several tributaries upstream of the intersection of the Callaghan and Cheakamus valleys (see Fig. 1). These occurrences are indicative of ephemeral higher stands of lava in the Callaghan Valley occurring at some point in the eruption and relate to a second feature. Immediately down drainage of the 45 m section of lavas exposed at Alexander Falls, the paleo-Callaghan Valley narrows to ≤100 m (Fig. 10C), which represents the potential for a natural blockage or impediment to flow (Figs. 10C and 11B). If we assume channelized flow in this narrow restriction, with depths of 1–2 m and depth-to-width ratios of 5, the cross-sectional areas would require velocities between 8 and 30 m/s to balance, for example, an eruptive flux of ∼160 m³ s⁻¹. These velocities are in excess of the velocities (1–5 m³ s⁻¹) arising from the average gradient (2–4°) in this part of the Callaghan Valley slope (using Jeffrey’s equation; viscosity of ∼500 Pa s). The consequence is that the eruptive flux of magma would be greater than the rate at which lava can flow through the restriction, and this could create the back pressure required to account for the forcing of lava to intrude subsidiary drainages.

The Cheakamus Valley contains one other example where interfingering lobes of EBp and BFp lava are found ∼1.2 km upstream of the junction of the Callaghan and Cheakamus valleys and ∼30 m higher in elevation (Fig. 1). This marks the northern extent of lava within the Cheakamus Valley. The upstream flow of lava reflects a backup of lava within the Cheakamus drainage caused by a high flux of lava exiting the Callaghan Valley, coupled to a bottleneck in the downstream flow as the Cheakamus Valley filled with basalt.

We suggest that the Cheakamus lavas represent a single continuous eruption event that occurred during the initial onset of the Fraser glaciation (≤25 ka). Precise dating of Pleistocene and younger mafic volcanic rocks by radiometric means (i.e., ⁴⁰Ar/³⁹Ar) remains a challenge due to magnitude of uncertainties that are commonly large enough (5–2 ka) to prevent discrimination of individual eruptions. We obtained three high-quality ⁴⁰Ar/³⁹Ar age estimates for the Cheakamus basalt lavas spanning 23.9 ± 7.9 to 13.1 ± 7.2 ka. Taken individually, it would be impossible to resolve whether they represent a single age of eruption or multiple eruptions of different ages; the data are permissive of both hypotheses. However, our paleomagnetic data established that a common paleomagnetic direction is shared by all Cheakamus lavas regardless of their stratigraphic position or geographic location (Fig. 8). The paleomagnetic pole direction does not fix the age of the volcanism, but it does restrict the total time represented by the lavas to less than the average rate of the Earth’s paleosecular variation (10^2–3 years). Furthermore, paleomagnetic data from thick sequences of lavas show no discernible systematic drift between the oldest and youngest lavas (Fig. 8B), suggesting even less time (<10² years).

On that basis, we can justify using the weighted mean of the three ⁴⁰Ar/³⁹Ar ages (15.95 ± 3.95 ka) as representing the best estimate of absolute age of eruption. That 1σ-weighted mean age window of 19.9–12 ka (Fig. 9) falls within the last glacial maximum at ∼16 ka (Clague and Ward 2011). At that time, ice thickness in the Fraser Lowland (∼100 km south of our study area) was at a maximum (∼2 km) and most definitely fully occupied the Cheakamus and Callaghan valleys. Our field study of the Cheakamus lavas shows the subaerial nature of the eruption and provides compelling field evidence for emplacement of the Cheakamus lavas immediately prior to peak Fraser glaciation. The evidence is diverse and includes (i) lavas inundating and filling the Callaghan Creek and Cheakamus River systems, indicating largely ice-free conditions; (ii) pervasive glaciation of the surfaces of the youngest lavas; and (iii) widespread distribution of glacial drift overlying Cheakamus lavas. In Fig. 9, we have compared the probability distribution function for the weighted mean ⁴⁰Ar/³⁹Ar ages to glacial events recorded in southwest British Columbia. We suggest an optimal eruption age of between 20 and 18 ka, perhaps, coinciding with the Port Moody interstade (Lian et al. 2001). While there is evidence for ice-free conditions in the Fraser Lowland at this time (e.g., Hicock and Armstrong 1985; Hicock and Lian 1995), there are no previous studies suggesting ice retreated past the Cheakamus Valley. However, rapid rates of ice advance and retreat are documented in other studies of the field area during the Fraser glaciation (e.g., Friele and Clague 2002). Our preferred window and duration of eruption is at odds with previous interpretations that invoked multiple periods of eruption preceding and coinciding with Fraser glaciation (e.g., Mathews 1958). Our field, paleomagnetic, and radiometric data require different durations, timings, and mechanisms of eruption for emplacement of the Cheakamus lavas.

Paleoenvironmental information is preserved in the overall morphology of lavas and through the distribution, orientation, size, and density of cooling joints. The lavas that infill the Cheakamus and Callaghan valleys preserve features suggesting subaerial eruption and lava emplacement into three different lava–paleoenvironment interactions, including ice-free conditions, ice-contact conditions, and water or meltwater infiltration.

The eruption of the Cheakamus basalts sourced from high elevation (∼1380 masl) within an ice-free Callaghan Valley. Localized restrictions in the Callaghan Valley created back pressures allowing lava egress into tributary drainages (Fig. 11B); this caused long-term damming of the paleo-Callaghan drainage to form the present-day Conflict and Callaghan lakes and diversion of these subsidiary drainages (Fig. 1). Indeed, the present Callaghan and Cheakamus Valley geomorphology is heavily controlled by the eruption of the Cheakamus basalts.

As the eruption continued, lavas entered the Cheakamus Valley, encountering the paleo-Cheakamus River. The lava flux was high enough to overwhelm the river system (Figs. 11B and 11C), and most lavas form subaerial, valley-wide, sheet flows (i.e., broad, laterally extensive blankets of lava) with well-developed decimetre-scale, lower colonnades capped by thick entablatures. Upper colonnades are almost always lacking, which could indicate that an upper colonnade did not form or was very thin and removed by glacial erosion. Notably, lower stratigraphic lavas exposed in the Brandywine Falls section were not glaciated (Fig. 2D) and also lack an upper colonnade.

The lavas within the Cheakamus Valley extend ∼10 km downstream and were likely fed by lava tubes. Partial collapse of these tubes has created a chaotic distribution of small lakes and potholes characterizing the surfaces of lavas within the Cheakamus River valley and the upper reaches of the Callaghan Valley. The emplacement of the basalt lavas into the Cheakamus Valley altered the path and incision rates of the paleo-Cheakamus River (e.g., Howard et al. 1982; Ely et al. 2012), which is today entrenched at the lava–bedrock contact on the eastern valley wall.

As the basalt lavas entered and progressively filled the Cheakamus Valley, they blocked the drainage system. The attendant consequence of this damming of the paleo-Cheakamus River was to build an ephemeral standing body of water upstream of the lava accumulation (e.g., Ely et al. 2012; Figs. 11C and 12). The northern extent of lavas within the Cheakamus Valley, described above (i.e., lavas north of Cheakamus–Callaghan junction), comprises complex, interfingering lobes of lava (Fig. 2A). These lavas mark the most upstream extent of lava and are remnants of the lava that dammed the paleo-river. These outcrops inform on the dynamic interplay between lavas inundating the Cheakamus River valley, damming the river, and creation of a syn-eruption upstream lake (Fig. 12).

In several places, outcrops comprise tabular sheets and lobes of lava (i.e., EBp vs. BFp) separated by fluvial to glacio-fluvial deposits (Fig. 2A). Commonly, the base of the overlying lava (i.e., BFp lava) is partially palagonitized (Fig. 2A) and features large (∼10+ cm) ragged vesicles and spiracles, suggesting lavas were emplaced on a wet, unconsolidated sediment-laden surface (i.e., Fuller 1931). In several instances, the lava has squeezed-up underlying sediment forming pronounced irregularities along the lava–sediment contact. The palagonitization and soft sediment deformation are important evidence for the short interval of time between sediment deposition and emplacement of subsequent lava(s).

The discontinuous lenses of sediments used to separate EBp from BFp lavas inform directly on the syn-eruptive competition between lava filling the Cheakamus Valley and the rise of the upstream lake level. The sediments are an indicator that at one point the lake level rose sufficiently to overtop the lava dam, creating a transitory sediment-bearing fluvial system. For as long as the lake level was maintained, it carried and deposited sediment across the lowest levels of the lavas spanning the Cheakamus Valley. None of the lava sequences in the Callaghan Valley host such interflow sediments, indicating that the damming-and-overtopping dynamic was limited to the Cheakamus Valley drainage system and did not occur in the Callaghan Valley.

The temporary flooding across the surfaces of the valley-filling lavas also explains some of the macroscopic textural properties of the EBp lavas in the Cheakamus Valley. The EBp lavas have sheet-like geometries and, although they have well-developed lower colonnades consistent with a mainly subaerial origin, they commonly lack an upper colonnade and always have thick entablatures featuring finer-scale and chaotic to fanning columns (Figs. 2B and 2D). These features reflect accelerated cooling of the lava and are commonly ascribed to ingress of coolant in the form of water or steam (e.g., Long and Wood 1986; Lyle 2000; Forbes et al. 2014). The pervasive distribution of water-infiltration features in the lavas suggests multiple episodes of localized damming of the river by lava and subsequent overtopping of the dam and flooding of the lava surfaces up to 11 km downstream (Figs. 11C and 11D). What is clear is that the processes of lava inundation, damming, and flooding are strongly coupled in time because the water infiltration occurs while the subaerially-emplaced lavas are still cooling.

Given the age established for the Cheakamus lavas, it is reasonable to consider to what extent ice may have occupied these drainage systems during the eruption. There is no evidence to suggest that ice occupied the Callaghan Valley. In fact, the distribution of lava and their characteristics (i.e., no glaciovolcanic features/lithofacies) would strongly suggest the Callaghan Valley was ice-free throughout the eruption. Similarly, the paleo-Cheakamus River valley must have been mainly ice-free. Most of the valley is filled wall to wall by flat lying, sheets of subaerial basalt lava, and there is a complete lack of glaciovolcanic features and lithofacies (e.g., Lescinsky and Fink 2000; Edwards and Russell 2002; Smellie and Edwards 2016).

This is true everywhere except at the Northern Terminus of lavas within the Cheakamus Valley (Figs. 1 and 11B–11D). These outcrops clearly record damming of the Cheakamus River and syn-eruptive overtopping and flooding of the lava dam (described above). However, there are also features that are not fully explained by those events. First, at this location within the Cheakamus Valley, the basalt lavas are found only on the west side, and even though they are at a high elevation relative to the valley floor, there are no equivalent outcrops of basalt lava on the east valley wall (Fig. 12). The implication is that lava flowed up the drainage system over 1 km and to elevations higher than the valley floor but could not cross the valley to reach the opposite wall. We suggest that ice partially occupied this part of the Cheakamus River valley but not further downstream. The corollary to this is that the Cheakamus River was glaciofluvial in origin, which is also consistent with the nature of some of the sediment interbedded between EBp and BFp lavas (Fig. 2A).

One other line of evidence derives from the lavas themselves. Large, overthickened (∼20 m), high-standing, vertical masses of lava exposed in the west bank of the present-day Cheakamus River exhibit highly exotic columnar jointing, including cube jointing, fanning to radially oriented columns, and, especially, horizontally oriented fine-scale (∼10 cm) columns that point out towards the valley center (e.g., Lescinsky and Fink 2000; Edwards and Russell 2002; Smellie and Edwards 2016). These steep vertical faces of lava record accelerated cooling and highly transient cooling histories and could be the result of lavas building up against localized ice masses and infiltration of meltwater (e.g., Mathews 1958; Edwards and Russell 2002; Smellie and Edwards 2016; Hodgetts et al. 2021). This would suggest that ice was present at or near the Callaghan–Cheakamus junction and worked in conjunction with the lavas to block and dam the drainage (Figs. 11C, 11D, and 12). In this case, much of the water accumulated in the ephemeral lake would be meltwater from glaciovolcanic interaction.

The end of the eruption is marked by the Upper Callaghan phase lavas (Fig. 1) that were emplaced subaerially and travelled at least 5 km down ice-free drainages (Fig. 11E). Subsequently, glaciers coalesced and advanced down the Callaghan Valley and into the Cheakamus Valley. The glacial advance eroded the vent area producing a flat, striated plateau, locally stripping the primary surfaces of exposed UCp, BFp, and EBp lavas and depositing 1–3 m of till. The till is still present within the Callaghan Valley but is mainly absent at lower elevations beginning at the junction of the Callaghan–Cheakamus valleys and continuing downstream.

The abrupt absence of till within the Cheakamus Valley (Blaise-Stevens 2008) indicates removal by a post-Fraser glaciation erosional event (i.e., 11 ka: Clague and Ward 2011). One possible explanation is erosion by a major glacial outburst flood event (e.g., Clague et al. 2020) that is suggested by the presence of (i) well-developed, smooth-walled drainage channels and escarpments up to 25 m deep incised into the basalt (Blaise-Stevens 2008), (ii) high-standing, scab-like remnants of BFp lavas (i.e., basalt “eskers”; Mathews 1948, 1958), (iii) irregular sculpted lava surfaces (Russell et al. 2007), and (iv) the misfit stream, amphitheater-headed canyon at Brandywine Falls (Fig. 2; Russell et al. 2007) and streamlined, butte-like morphology of Daisy Lake Island. This would suggest that the present-day morphology and distribution of the basalt lavas within the Cheakamus Valley results, in part, from a major flood event during waning of the Fraser glaciation. In this scenario, the landforms that Mathews (1958) postulated to be “basalt eskers” may be remnant scarps from more extensive sheets of lava excavated by one or more glacial outburst floods.

The Cheakamus basalt lavas sourced from a vent situated near the present-day Conflict Lake during the advance of the Fraser glaciation. Minimal petrographic and geochemical diversities suggest that the basalts were erupted from a common batch of magma that underwent minor chemical differentiation involving Ol and Pl phenocrysts and entrainment of xenocrystic or antecrystic Pl. Our paleomagnetic studies show all lavas to share a common paleomagnetic direction (346.18°/72.97°; α₉₅ = 1.3°), indicating that the eruption occurred within a single paleomagnetic moment (<2000 years) rather than over a protracted time as suggested by previous studies. Furthermore, a close examination of paleomagnetic results from the base and top of thick stratigraphic sections shows no systematic drift in direction, suggesting that the duration may be even shorter. The weighted mean of three ⁴⁰Ar/³⁹Ar dates provides an absolute eruption age of 15.95 ± 3.95 ka (1σ), which combined with field evidence (i.e., glaciated surfaces) is consistent with the pre- or early-stage Fraser glaciation, possibly during the Port Moody interstade (i.e., ∼20–18 ka).

The Cheakamus basalts were erupted subaerially into major mountainous drainages typical of the Coast Mountains. The lavas advanced 14 km from a high elevation source (1380 masl) down an ice-free Callaghan Valley to enter a mainly ice-free, paleo-Cheakamus River valley at 506 masl. The lavas entering the Cheakamus River valley backed up and flowed upstream ∼1 km from the valley junctions and flowed downstream for another 12 km, causing damming of the drainage system. The lava dam produced a substantial (≤30 m deep) ephemeral, upstream standing body of water that periodically overtopped the lavas to flood their surface. The overtopping and flooding are manifest by discontinuous lenses of fluvial–glaciofluvial sediments and thick entablatures within the lavas indicative of accelerated cooling driven by syn-eruptive (i.e., cooling) water infiltration. The distribution and properties of lavas occupying the northern end of the Cheakamus Valley suggest that localized masses of ice partially filled this part of the drainage and contributed to the damming of the paleo-river system. Evidence for a major post-eruption and post-Fraser erosional event may indicate the occurrence of one or more glacial outburst floods partially accounting for the present-day distribution and character of Cheakamus basalts.

Author contributions: Geological mapping and sampling was carried out by AB. RWB assisted with the paleomagnetic sampling program and carried out the paleomagnetic measurements. The project was originally conceived by JKR, and all authors were involved in the interpretation and implications of the datasets, as well as the generation of the figures and the writing of the paper. Several high-resolution field photographs were kindly provided by Steve Quane. The manuscript benefited immensely from critical reviews provided by Dave McGarvie and an anonymous referee.

Data generated or analyzed during this study are provided in full within the published article and its appendices.

Conceptualization: JKR

Data curation: AB

Formal analysis: AB, JKR

Funding acquisition: JKR

Investigation: AB

Methodology: AB, JKR, RWB

Resources: JKR, RWB

Software: RWB

Supervision: JKR

Visualization: AB

Writing – original draft: AB, JKR

Writing – review & editing: AB, JKR, RWB

This research was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery grants to JKR and RWB and NSERC PGS-D awarded to AB.

The Cheakamus basalts show two varieties in terms of their alternating field and thermal demagnetization properties. Samples 212 and 315 are consistent with most samples revealing magnetizations residing predominantly in fine-grained single-domain magnetite crystals (i.e., “hard magnetization”). The thermal demagnetization curves (left) for these samples have the classic “shoulder” shape where intensity does not change much with stepwise demagnetization until ∼500–580 °C. Beyond that heating step the magnetization drops rapidly, typical of fine-grained single domain magnetite. In contrast, samples 017 and 131 show a mix of fine-grained single domain magnetite and coarser-grained multi-domain magnetite (i.e., “soft magnetization”). In these samples, magnetization drops more or less uniformly with stepwise demagnetization that is typical of a predominance of coarser-grained multi-domain magnetite.

Three holocrystalline, non-vesicular samples of the Cheakamus basalts were collected for ⁴⁰Ar/³⁹Ar analysis and age determinations (Table D1). Two samples were submitted to the Argon Geochronology Laboratory at Oregon State University. The other sample was submitted to WiscAr Geochronology Labs at the University of Wisconsin–Madison. Samples were washed and hand-picked for phenocryst-free whole rock chips of crystalline groundmass. The monitor mineral TCR‐2 with an age of 28.619 Ma (Renne et al. 2010) was used to monitor neutron flux and calculate the irradiation parameter (J) for all samples. Instrument detectors simultaneously measured ⁴⁰Ar/³⁹Ar, ³⁶Ar/³⁹Ar, and ³⁷Ar/³⁹Ar ratios through cycles of beam switching. Ages were calculated relative to the 1.1864 Ma Alder Creek sanidine standard (Jicha et al. 2016). Argon decay constants and isotopic abundances are after Min et al. (2000).

Results are summarized in Table D1 and Fig. D1. The plateau ages are preferred over the inverse isochron ages due to higher precision for sample AW-18-22. Samples AB-20-33 and AB-20-45 have similar plateau and inverse isochron ages. The spectrum provides an acceptable plateau age if three or more consecutive gas fractions represent at least 50% of the total gas release and are within two standard deviations of each other (mean square-weighted deviation less than 2.5). Spectrum plateau ages are quoted to 1 s (66.7% confidence) and calculated using the constants of Renne et al. (2010). The integrated age is the age given by the total gas measured and is equivalent to a potassium–argon (K–Ar) age. One sample (AB-20-33) had an integrated age within error of the plateau and isochron ages (Table D1).