Where ice gave way to fire: deglacial volcanic activity at the edge of the Coast Mountains in Milbanke Sound, BC Open Access

Glaciovolcanism occurs where local glaciers or regional ice caps are present during volcanic activity, with numerous examples in British Columbia (BC) (Mathews 1947; Grove 1974; Hickson 2000; Edwards et al. 2002; Hickson and Vigouroux 2014; Wilson and Russell 2020) as well as Antarctica, Iceland, Chile, Patagonia, New Zealand, and Kamchatka (e.g., Vinogradov 1981; Bösken and Nett 2013; Conway et al. 2015, 2023). Russell et al. (2014) explained volcanic facies, processes, and genetic models for subglacial volcanism and tuya formation. There are compelling and complex questions concerning causal links between glaciation, ice loads, mantle flow, melt generation in the asthenosphere, and melt delivery to the surface (Hall 1982; MacLennan et al. 2002; Schmidt et al. 2013; Watt et al. 2013; Rawson et al. 2016; Mora and Tassara 2019; Wilson and Russell 2020).

The thickness of the Cordilleran Ice Sheet (CIS) over western Canada and Alaska at the last glacial maximum (LGM, ∼21 cal ka BP) of >1 km (Dyke et al. 2003; Dyke 2004) is tantamount to loading and removing more than 300 m of rock within a few thousand years (>10 cm/year), at least an order of magnitude higher than current tectonic strain far from the continental margin (Brillon et al. 2018). Deglaciation enhances rates of mantle decompression melting and volcanism, as interpreted for Iceland (MacLennan et al. 2002; Schmidt et al. 2013) and Southeast Alaska (Praetorius et al. 2016; Supplementary File 1). Basaltic magmas form by decompression melting of mantle peridotite (e.g., Hirose and Kushiro 1993). Sufficient glacial loading and unloading might move enough asthenosphere to generate small volume magmas for this within-plate setting. The glacial rejuvenation of faulting is also likely of great importance to magma transport and the siting for deglacial volcanoes.

This study concerns a 43 km line of small Quaternary volcanic piles exposed around Milbanke Sound on Canada’s central Pacific coast. They lie in the unceded territory of the Klemtu-Kitasoo/Xaixais and the Bella Bella-Heiltsuk peoples. Kitasu Hill (KH) and four other eroded, small volume basaltic scoria piles and flows on nearby islands were known from regional surveys, and variously called the MSV (Souther 1990), Lady Island Formation (Dolmage 1922, 1924), and Lake Island Formation (Nelson et al. 2011b; 2014). The volcanoes were thought to be recent and incompletely mapped or investigated by previous workers due to remote access and dense rainforest.

We sampled two of the MSV (Figs. 1a–1e). Volcanic rocks and volcano-glacial sediments of KH on Swindle Island were sampled to ground truth a new LIDAR (LIght Detection And Ranging) survey of KH and to characterize its petrology (Bednarski and Hamilton 2019). Outcrop was very restricted to the beach, a densely forested ravine and the occasional root zones of fallen trees. We discovered McGregor Cone (MC) during routine marine surveying from recent multibeam bathymetry. It lies atop Finlayson Moraine beneath Milbanke Sound. These volcanoes sit on glacially eroded surfaces and contain coeval glacial and volcanic materials. MC has foraminifera that lived on the seabed at the time of deposition. We interpret that both volcanoes formed during deglaciation and offer details about the sequence of events and their controlling processes.

This study is not a typical case history of two small volcanoes; rather it uses aspects of the volcanic geology to infer the glacial and relative sea level (RSL) history. We used a multidisciplinary approach to integrate the volcanology, geochemistry, geochronology, and geomorphology from a LIDAR survey of Swindle and Price Islands, a new digital elevation model (DEM) and four marine mapping and sampling cruises in Milbanke Sound. We examine volcano-glacial and postglacial facies, glacial and eruptive history, and the type of volcanic products in a sequence stratigraphic framework. We provide an RSL curve from the direct dating of benthic foraminifera to constrain the sequence of eruptive events and sedimentary facies during the waning stages of the last ice age. Our study provides details about enigmatic small-volume intraplate volcanism in western Canada and offers clues into connections between deglaciation and volcanism in such settings. We propose a model to explain the relevance of these small deglacial volcanoes to rest of the Cordilleran region.

The geomorphic interpretation presented here relies on high resolution topographic and bathymetric data. Analysis of the KH volcano was by LIDAR imagery (provided by Brian Menounos of UNBC, technical and processing details in Bednarski and Hamilton 2019). The DEM for Canada’s Pacific coast (Kung et al. 2023) includes bathymetry from 2010 that remained uninterpreted until 2020. The multibeam data have a vertical resolution of better than 2 cm and a slope resolution better than 0.25°.

The new bathymetry in the DEM revealed the steep symmetric shape of MC (Figs. 1d, 2b, and 3). MC sits in a gap in the line of previously documented volcanoes along the Principe Laredo Fault. It rests on a recessional moraine across the mouth of Finlayson Channel. Subsequent campaigns (CCGS Vector cruises 2015003PGC, 2020OSD2020, 2021003PGC, and 2021006PGC) used CHIRP (Compressed High-Intensity Radiated Pulse) subbottom seismic profiles, bottom photographic traverses, remotely operated vehicle (ROV) videos, grabs of seabed volcanic rock, glaciomarine sediment and postglacial mud, and gravity cores of glaciomarine and modern sediment at diverse sites to confirm MC’s volcanic origin and to study its composition, stratigraphy, and timing.

Whole rock samples (5–30 g) from KH and MC were trimmed by saw to remove alteration, washed, dried at 80 °C, and submitted to commercial labs for pulping, Loss-on-ignition and major and trace element analysis by inductively coupled plasma mass spectrometry (ICPMS) or emission spectroscopy (ICPES) with replicates. The KH samples were analysed by Acme Labs, Vancouver (Bureau Veritas Commodities Canada Ltd.; Bednarski and Hamilton 2019). The MC samples were analysed by ALS Global in Vancouver using their procedures: BAT-01, PREP-31, and CCP-PKG 01. Analyses included a lithium borate fusion ICPMS and ICPES/MS with a modified aqua regia digestion (1:1:1 HNO₃:HCl:H₂O). The ALS analyses included replicates of two KH samples for interlab comparison.

We also geochemically analysed samples of MC tephra from washed and dried, sand-sized separates (>63 µm) at two levels in two cores. Proximal core 93B was taken downslope from the cone’s amphitheatre at 131 m below sea level (mbsl). Distal core 124 was 1.5 km away, not directly downslope to avoid possible tractive transport and upsection reworking of shards. Clean, glass fraction concentrates were hand separated under a binocular microscope to obtain 8–85 mg. At the University of Alberta Tephrochronology Lab, glassy particles were further hand-picked and bulk mounted in acyclic pucks with epoxy, and polished and carbon-coated for electron probe microanalysis. Points were picked under backscatter imaging to avoid microlites and other mineral contamination. Analyses for major elements were carried out on a JEOL 8900 SuperProbe by wavelength-dispersive spectrometry using a 15 KeV voltage, 6 nA current, and 5 µm beam. A time-dependent intensity correction was used to correct for potential Na-loss caused by the highly focused beam (Donovan et al. 2015). Two secondary glass standards, Lipari obsidian ID 3506 (Kuehn et al. 2011) and BHVO-2G (Jochum et al. 2005), were run concurrently to track analytical accuracy and calibration stability. Individual analyses normalized to 100% are presented as weight percent oxides along with averages and standard deviations (Table S5).

MC airfall tephra were directly dated by wet sieving and hand-picking benthic foraminifera, which lived on the seabed at the same time the tephra and glaciomarine diamicton fell from floating ice. These dates fix the timing of the final glaciomarine strata and waning volcanic activity in the MSV.

Samples from sediment cores were wet‐sieved through a 63 µm screen, and 4–12 mg of foraminiferal specimens from each sample were identified and handpicked for accelerator mass spectrometry (AMS) measurements at the National Ocean Sciences AMS facility at Woods Hole Oceanographic Institution. Radiocarbon ages were calculated using the method of Stuiver and Polach (1977).

Calibrations used OxCal Version 4.4 (Bronk Ramsey 2009) with the Marine20 calibration curve (Heaton et al. 2020; 2022). Marine20 is explicitly valid only for surface waters (<100 m), but only benthic foraminifera, which typically require larger carbon reservoir corrections, were identified in our samples. Care was taken to avoid including infaunal taxa, which grew with unknowably old carbon. The current water depths of our samples are 131–274 mbsl; however, sea level was much lower at the time of deposition. Thus, we applied the regional ΔR₂₀ corrections typically applied to planktic foraminifera, using the Northwest Coast of North America averages of Schmuck et al. (2021). The Holocene-age samples from Station 123 use a ΔR₂₀ correction of 145 ± 165 years, while the rest use the Bølling–Allerød interstadial age (14–10.7 cal ka BP) correction of 575 ± 165 years.

The MSV straddle the Principe Laredo Fault (Fig. 1a), a right lateral Miocene crustally penetrative shear zone along the northeast side of Queen Charlotte Basin (QCB; Rohr and Dietrich 1992). MSV lie at the Insular Belt’s inner edge (Supplementary File 1) and was its last part to deglaciate, based on a succession of progressively deeper and older drowned moraines to the southwest (Fig. 2a). Milbanke Sound has four distinct physiographic subregions (Figs. 1a, 2b, and 3), including (1) the Milbanke Strandflat (Slaymaker and Kovanen 2016), (2) flat-floored shelf-crossing troughs of Laredo and Milbanke Sound, (3) dissected alpine uplands of the outer fjord ranges, and (4) overdeepened fjord channels. The erosional morphologies above (75 masl) and below sea level (200 mbsl), along with volcanic and sedimentary facies, reveal the geological history since the Wisconsinan ice age’s end in Milbanke Sound.

The Milbanke Strandflat, spanning southwest Swindle Island through the southwest Don Peninsula, contains all the subaerial MSV (beige interval 0–75 masl in Fig. 2b). It is a flat, poorly drained rocky terrain with expanses of muskeg, lying mostly below 30 masl, in contrast to the rugged fjord landscape immediately to the northeast. Subglacial erosion must have been a significant factor because exposed bedrock surfaces show intense scour and distinct lineaments; near KH on Swindle Island and Price Island, the glacial incisions trend southwesterly (fig. 12 of Bednarski and Hamilton 2019). The strandflat exposes and truncates diverse bedrock types (Baer 1973), especially Early Cretaceous plutons and thin metasedimentary infaulted blocks of the Devonian Matheison Channel Formation (Fig. 1a) of the Paleozoic Alexander Terrane (Nelson et al. 2011a, 2011b, 2014). It contains thin remnant flows of Masset/Bella Bella Formation (Nelson et al. 2014) of previously much thicker (>5 km) rift basin fill. Rohr and Dietrich (1992) show postrifting (Neogene) movement here along the Principe Laredo Fault. This strandflat is an exhumed Tertiary unconformity and a compound Quaternary erosional surface.

Swindle Island north of Higgins Passage has a fluvioglacial outwash delta (fig. 12 of Bednarski and Hamilton 2019), built from ice to the NE when the marine limit was at ∼65 masl. If the 75 masl contour approximates a former sea level after deglaciation but prior to rebound, the only exposed land areas would have been the grey and tan areas (Fig. 2b) including KH, an archipelago of small islands on eastern Price Island, and the outer fjord ranges in the NE quadrant of the map. Queen Charlotte Sound (QCS) extended across the strandflat (Fig. 2b, tan and brown) and eroded it to wave base with the fetch of the open Pacific Ocean.

The low-lying (∼270 mbsl) troughs of Milbanke Sound and Laredo Sound comprise the second young physiography. The north-trending linear coastlines of Price Island and adjacent islands (Figs. 1a and 2) suggest Milbanke and Laredo sounds are underlain by grabens. These local flat-floored sounds were incised across the Milbanke Strand Flat during one or more glaciations. Unlike fjords, they do not have U-shaped profiles and seem to have been partially formed by flowing water with its tractive base level beneath or beyond the ice. These shelf-crossing troughs were partly infilled by lateral and recessional glacial moraines followed by distal glaciomarine deposition and finally more restricted drifts of postglacial hemipelagic mud.

A third physiography is the higher elevation (100–800 masl) rugged alpine peaks and ridges of the outer fjord ranges on Swindle, Dowager, Susan, and Roderick Islands and the northern portion of the Don peninsula. Unlike the smoothed whaleback profiles from the CIS, this more rugged local topography has horns, aretes, and cirques that are formed by alpine periglacial processes as supported by Holocene alpine glacial records to the north and south of this region since 8.9 ka (Mood and Smith 2015; Darvill et al. 2022). This sculpting incised down to ∼100 masl in local U-shaped glacial valleys.

The fourth, deepest, and most restricted physiography encompasses the overdeepened channels incised into the seabeds of the fjords and sounds, shown in dark blue and purple (Fig. 2b, >275 mbsl). These regions include Laredo Sound west of Price and Swindle Islands, Finlayson Channel, Mathieson Channel, and Spiller Channel (all three are deeper than 650 mbsl). The deeps lie up-fjord from ice grounded on recessional moraines, with floating ice shelf conditions behind them. Hydraulic erosion during outburst floods generated these overdeepened basins like the one at Laskeek Bank (Shaw et al. 2020) and elsewhere at the mouths of Coast Mountains fjord such as Douglas Channel circa 11.9 cal ka BP (Shaw et al. 2017). The overdeepening of fjords, with hydraulic fluting in the downflow direction, is common across the continental shelf of western Canada, when ice dams floated off to release built up water and ice from the upper fjord system, like the glacial outburst floods of 2016 in Jakobshavn Fjord, West Greenland. The local hydraulic fluting scored across moraines, and eroded into weaker faulted and folded bedrock, beneath Laredo Sound, Principe Laredo Sound, and Finlayson Moraine at the southern end of Finlayson Channel (Fig. 3). The seabed of these channels is mostly bare rock or has a thin cover of glaciomarine silt and less than a few metres of Holocene mud. Ice withdrawal, glaciomarine sedimentation, and deglacial failures in the Douglas Channel–Kitimat fjord system just to the north (Bornhold 1983) began about 13 cal ka BP. There are few modern sediment sources to deliver fine sediments to Milbanke Sound during the Holocene transgression and highstand conditions since ∼7 cal ka BP.

Tectonics have influenced glaciation as much as the climate has along the Northeast Pacific Ocean (Mann and Hamilton 1995). The Queen Charlotte–Fairweather Fault system’s transpression since Miocene accompanied the uplift of the Alexander Archipelago in Southeast Alaska (Plafker et al. 1976) and the BC Coast Mountains. The western lobe of the CIS arose multiple times since late Miocene (Eyles et al. 1991; Lagoe et al. 1993) and attained maximum thickness in these high elevation coastal regions. Here, massive outlet glaciers (Goldthwaite 1987) formed deep troughs, which cross the continental shelf and incise well below current sea level, through the Alexander archipelago. A series of recessional till (moraine) tongues (Josenhans et al. 1993) are below (265–400 mbsl) eustatic sea levels along Moresby Trough (Fig. 2a). Drowned lacustrine deposits (Josenhans et al. 1997) and soils containing tree roots lie beneath QCS (Barrie et al. 1993). Complex glacial facies, bedforms, and erosional morphologies indicate a dynamic deglacial and tectonic history, requiring both substantial rebound and a glacial forebulge. The variability in RSL curves across the margin (Shugar et al. 2014) suggests that local faults control relative timing for early rapid rebound and later more gradual forebulge collapse.

In contrast to the Laurentide Ice cap of Eastern North America, which sat on a thick, cold, ancient craton (Peltier 1998, 2004), this region of the Coast Mountains and Insular Belt was a rift basin throughout middle Tertiary and still has thin crust, shallow, warm, low-viscosity asthenosphere whether due to the recent rifting or to the overriding of warm recently upwelled Pacific mantle (Lewis et al. 1991, 1997; Lowe and Dehler 1995; James et al. 2000; Hyndman 2010, 2017; Hyndman and Canil 2021) and abundant throughgoing faults (Rohr and Dietrich 1992). This active tectonic setting exerted a strong influence on the nature of glacio-isostatic motions both locally, with faults to accommodate rebound on a small scale, and regionally through unusually active, warm, and shallow asthenospheric mantle flow.

At the Wisconsinan onset, prior to local LGM, increasing ice load in the Coast Mountains sent valley glaciers down the fjords to coalesce as piedmont glaciers in Hecate Strait and QCS, with a separate lobe traversing Dixon Entrance (Barrie and Conway 1999). As the thickness and weight of ice built up locally across the shelf, local glacial isostatic depression partly countered the forebulge from the growing CIS. RSL controlled the CIS Pacific limit (Mann and Hamilton 1995; Darvill et al. 2022). LGM and ice retreat were out of phase and diachronous along the coast, wherein the CIS reached its LGM and initiated deglaciation earlier in Alaska than coastal BC. Mann and Hamilton (1995) put LGM for Alaska’s peninsular tidewater glaciers between 23 and 16 ka ¹⁴C BP (28–19 cal ka BP) but somewhat later, 15–14 ka ¹⁴C BP (18–17 cal ka BP), farther south along the BC coast. More recent ¹⁰Be dates on erratics and bedrock exposures (Darvill et al. 2022) interpret ice retreat to the inner coast (including Milbanke Sound) by 18–16 ka, with thinning to 13 ka as ice withdrew from marine areas and became fully terrestrial. Barrie and Conway’s (1999) analysis of shelf sediments from Dixon Entrance, Hecate Strait, and QCS also support this timing. Dates on deglacial deposits of the continental margin of Vancouver Island, south to the Strait of Juan de Fuca, show deglaciation by 14 cal ka BP (Hamilton et al. 2015).

Recessional moraines across the shelf (Fig. 2a) and in Milbanke Sound grounded near local sea level. The moraine toes and the base levels they imply are much deeper than ∼130 mbsl, the LGM far-field eustatic minimum (Lambeck et al. 2014). We infer that moraines were uplifted to equal global eustatic sea level at the time of grounding. Josenhans et al. (1993) called this distal uplift away from the centre of the ice sheet the “forebulge” and used it to explain deep till tongues and drowned lacustrine deposits in Hecate Strait. Massive ice on the Coast Mountains displaced asthenosphere to the west (against regional tectonic flow) and uplifted a flap of thinned crust across the Insular Belt as a forebulge (Josenhans et al. 1993). QCB’s crust is thin, weak, and extensively faulted (Rohr and Dietrich 1992), ruling out crustal flexure.

MSV’s five small basaltic centres rest on glaciated bedrock and are well preserved, suggesting they were postglacial and (or) Holocene (Baer 1973). We note that the MSV sits on the southeast tip of the Alexander terrane between the Principe Laredo Fault and Bonilla Island Fault (Rohr and Dietrich 1992). On our LIDAR DEM spanning Milbanke Sound (Figs. 1a, 2b, and 3), subglacial erosion has enlarged contacts, faults, and joints cutting Cretaceous plutons (figs. 3 and 12 of Bednarski and Hamilton 2019). Using these data, we tentatively relocate the southeast end of the Principe Laredo Fault as curving to the east across Dowager Island (Fig. 1a). The linearity of Mathieson Channel and Spiller channel on the DEM, and mismatch of bedrock plutons across them, suggest north northeast-trending faults control those channels (Nelson et al. 2014). Helmet Peak (Fig. 1e) lies exactly where Mathieson Channel intersects the projected Principe Laredo Fault. The MSV trend probably reflects a deeper mantle source region or set of deep crustal magma chambers beneath this faulted uplifted plutonic crust.

Tertiary reactivation of crustal faults occurred as the Queen Charlotte rift basin developed. (Fig. 1a; Rohr and Dietrich 1992). The Principe Laredo Fault zone, across the heads of Milbanke and Laredo Sounds, is part of a transtensional system that was formed during middle Tertiary dextral extensional faulting. This resembles the Tertiary Masset Formation filled grabens beneath Hecate Strait and QCS. The Principe Laredo Fault has 80 km of dextral offset since mid-Cretaceous, based on displacement of the northern end of the Banks Island plutonic suite (Nelson et al. 2014). Irrespective of any modern tectonic strain on these faults, as they cut the crust and thinned mantle lithosphere, they likely served as conduits for magma delivery. Given the large- and short-term glacial loading and unloading since LGM, many of these deep faults were likely reactivated, as ice loads changed as per neoglacial faulting in Puget Sound (Thorson 1996, 2000). These major faults allowed rebound to be accommodated locally, one fault block at a time (elevator style), rather than through some broader scale crustal flexure.

There is high heat flow for this region, greater than 70 mW/m² (Lewis et al. 1991, 1997) due to thinned lithosphere and a shallow warm upper mantle (Lowe and Dehler 1995; Hyndman 2010, 2017; Hyndman and Canil 2021). Evidence for vertical fault motion since Miocene rifting includes erosion down to the basal nonconformity for QCB where the Masset volcanics (Bella Bella formation) overlie Mesozoic plutons, from Porcher Island and Aristazabal Island to the northwest (Woodsworth 1991) through to Dufferin Island in the southeast (Fig. 1a).

On land, four MSV volcanic centres with five outcrop areas span of 43 km (Fig. 1a). KH on Swindle Island has the most complex volcanic facies and is the best preserved, perhaps from a more prolonged history than is typical for monogenetic cinder cones. Baer (1973) described basaltic columnar flows on the east shore of Price Island 3 km south of the mouth of Higgins Passage. These flows are 6 m thick and dip 20° seawards (implying post-eruption tilt) and are likely the outcrops mapped by Dolmage (1922, 1924). Price Island Cone is rubbly and its flow near sea level is eroded (Souther 1990). The Price Island outcrops are not mapped on the Nelson et al. (2014) compilation; however, we tentatively locate the Price Island cone 3 km south of the Principe Laredo Fault (Fig. 1a) based on prominent resistant mounds expressed in the LIDAR DEM. With smoothing and possible wave base erosion up to ∼75 masl (Fig. 2b), it is likely that the Price was more severely eroded than higher elevation centres like KH (233 masl) and Helmet Peak (285 masl), which escaped the wave base erosion across the Milbanke Strandflat. Helmet Peak on Lake Island (Figs. 1a and 1e) is a prominent vent with pyroclastics, dykes, hyaloclastites, and large erratic plutonic blocks (Souther 1990), so it was likely emplaced subglacially as per other tuya-type features (Russell et al. 2014). From the DEM (Fig. 1e and fig. 13 of Bednarski and Hamilton 2019), its steep cliffs have been eroded either by late glacial outburst floods down Mathieson Channel or wave base erosion at higher RSLs up to 75 masl. Helmet Peak is erosionally separated from its distal pyroclastic breccias to the west on Lady Douglas Island and from a smaller unnamed island to the east, where any local sea level higher than 40 masl would suffice. The thin Dufferin Island flows (Souther 1990; Nelson et al. 2014) are low lying and eroded on the Milbanke Strandflat; their vent has not been identified.

Other late Quaternary Insular Belt volcanoes are scattered across the Alexander archipelago in southeast Alaska between the Mount Edgecumbe volcanic field (MEV) on Kruzof Island (Riehle et al. 1992; Praetorius et al. 2016; Grapenthin et al. 2022) and the edge of Dixon Entrance at the Canadian border. They include Tlevak Strait-Suemez Island Volcano near Prince of Wales Island (Brew 1990; Riehle et al. 1992), the drowned cinder cone near Misty Fjords National Monument (Karl et al. 2013; Kheiry 2013), and the submerged Addington basaltic volcanic field (Wilcox et al. 2019) in southeast Alaska. They all sit atop glaciated Alexander Terrane bedrock, directly against the Queen Charlotte Fault and the Pacific Plate. These late Wisconsinan lavas in Alaska, like MSV, are all proximal to uplifted and eroded Tertiary rift basin lavas associated with that penultimate and more voluminous volcanic episode for the Insular Belt, and all are outboard of the thickest CIS load on the Coast Mountains.

KH volcano (Fig. 1c, Supplementary File 1) is a nested cindercone complex with later backstepping flows on a squared-off steep-sided platform 2.3 km across and with a total volume of 0.32 km³. Bednarski and Hamilton (2019) recognized that KH started subglacially as the ice thinned. The geomorphic evolution of KH as revised here is presented in Fig. 4. KH initially built a shallow dipping set of volcano-glacial sediments with overlying pillows and pillow breccias to form a tuya, all within a subglacial or englacial lake. As ice withdrew, the lake drained and the edges of the tuya eroded forming squared off cliffs and beaches when local sea level was at 42 masl, early during postglacial rebound. The facies, samples (Table S1), and their geomorphic evolution interpreted from the new DEM, are presented in Supplementary File 2. Early deglacial conditions, subsequent erosion and falling sea levels during rebound make the geomorphic development in Milbanke Sound differ from distant unglaciated areas that only experienced the Holocene transgression.

Recent multibeam surveys revealed MC, a fifth vent in the MSV (Figs. 1d, 2b, and 3), atop Finlayson Moraine, at the southern end of Finlayson Channel. Milbanke Sound has a complex multistage erosional history. This composite surface was successively modified by glacial abrasion beneath the CIS, early higher sea levels, and stronger coastal storms after ice withdrawal but prior to rebound, deglacial outburst floods from farther up the fjords, a marine regression during postglacial rebound, a prolonged transgression since ∼15 ka and strong persistent tidal currents up to four knots. Sample locations for photos, grabs, and cores are listed in Tables S2a, S2b, and S2c. Further details are presented in Supplementary File 3.

We find that bedrock ridges along overdeepened Finlayson Channel have prominent peaks and ridges along their axes, whereas the slopes climbing Finlayson Moraine, along its deeper southeast edge, and the sides of MC are more rounded and separated by regions of pronounced hydraulic fluting (Fig. 5). Ice withdrawal from Milbanke Sound formed Higgins Moraine across the mouth of Higgins Passage, and Finlayson Moraine across the mouth of Finlayson Channel. On CHIRP lines, these steep-sided moraines have smooth profiles (Figs. 6a and 6c) and are acoustically chaotic, lacking internal reflections. Both moraine diamicton and bedrock show strong seabed returns and high intensity backscatter.

Finlayson Moraine extends across the channel with a well-defined 50 m high scarp along its southwest face, likely where ice grounded. The deeper northeast slopes of Finlayson Moraine are rugged and eroded by subsequent ice and hydraulic action. There are a few blocks of dense diamicton in front of Finlayson Moraine (Figs. 2b, 3, and 5) as if from mass flow events. The moraine shallows to 150 mbsl against Swindle Island and Pidwell Reef in Higgins Passage and against Dowager Island to the southeast, but there is a steep drop into Finlayson Channel descending to greater than 650 mbsl (Figs. 2b and 3). The upstream (northeast) slope of the moraine has chevron-shaped ridges, with limbs up to 500 m long, pointing southwest downflow. The southernmost chevrons cut across most of the moraine, with an iceberg-pitted and scoured surface at 170–120 mbsl and elongate hydraulic flutes, like those from glacier outburst floods. This erosion occurred when floating ice in Finlayson Channel was dammed against the terminal moraine but before MC formed.

Both moraines are built of diamicton and are onlapped by glaciomarine sediment, deeper than 130 mbsl. These massive to partially laminated silty beds contain thin sandy and fining upwards sequences, dense water-lain, organic-poor deposits with faceted grit, sparse marine fauna, and glacial dropstones at the seabed as per bottom photos (Figs. 7e and 7f). We interpret the glaciomarine facies to have been deposited under floating ice shelf conditions, beyond the calving receding ice front, by combined hemipelagic settling and tractive basal deposition from icebergs. Glaciomarine sedimentation rates near MC (Fig. 1d, Table 1) are 0.98–1.97 mm/year, like the distal deglacial sedimentation rates on the deep continental slope (Hamilton et al. 2015) or the foredeeps of the Fraser River (Hart et al. 1995). Onlapping glaciomarine sediment with airfall tephra onto MC’s lower portion demonstrates that the subaerial volcano formed earlier and continued to erupt explosively as it became drowned and draped by glaciomarine materials.

Strong CHIRP reflections are common where the seabed is diamicton (Figs. 6b, 6c, 7e, and 7f). Glaciomarine sediments can be distinguished from postglacial mud by strong and continuous, bottom hugging, internal reflections in seismic profiles (Figs. 6a and 6c). This glaciomarine unit usually has a hard lag surface, with strong backscatter, armoured by glacial ice-rafted debris (IRD), but it has internally stratified reflectors that drape underlying units. The moraines bottom photos show a lag of rounded cobbles and boulders, with a benthic crust. Grab samples have rounded and weathered coast plutonic rocks and occasional red oxidized volcanic bombs, blocks, and scoria on MC (Figs. 7a–7c).

The Higgins Moraine has a series of downward terraces from Price Island to the northwest. These formed as a regressive surface of marine erosion during a late Wisconsinan regression, due to local deglacial rebound and persistent forebulge conditions. Depths are greater than 180 mbsl on the Milbanke Sound side (southeast) shallowing to 130–150 mbsl atop the moraine and dropping to 260 mbsl on the upflow (northwest) side. This implies a minimum of ∼215 m of postglacial rebound on Higgins Moraine next to the Principe Laredo Fault and Swindle Island. Diamicton is exposed along its southeast side in Milbanke Sound, while most of the flat top of the moraine is thinly mantled by glaciomarine sediment and IRD (Figs. 7c, 7e, and 7f) from a subsequent transgression during floating ice shelf conditions, as local rebound deceased and the forebulge collapsed back to the northeast. The glaciomarine sediment mantling the Higgins moraine is iceberg-pitted with extensive ornate build-ups of glass sponge reefs, like other well-flushed, glaciated seabeds, now in hiatus up and down the west coast from Hecate Strait (Barrie and Conway 1999) to the Salish Sea.

Finlayson Moraine flattens out to the northwest end near 120 mbsl and is overlain by MC volcano, a broad conical feature (Figs. 3 and 5) with ∼26° side slopes and even steeper slopes (>40°) near the summit crater (Figs. 5 and 6b), like many subaerial cinder cones (Wood 1980). This cone is littered with oxidized scoria, blocks, bombs, and tephra (Figs. 7a–7c; Table S2A) whose size decreases downhill and away from the central vent in all directions. The volcano was subaerial and formed atop the Finlayson Moraine after it was hydraulically eroded and rebounded above local sea level.

The cone rises from ∼200 to 43 mbsl, with a basal diameter of 2.6 km and a volume of 1.1 km³. As per KH, the volume is small compared to monogenetic cinder cones worldwide (Wood 1980). The cone is breached on its southwest flank and has several irregular curvilinear ridges (Figs. 3 and 5). The two ridges that run directly downslope from the breached crater edges from about 50 to 125 m on the southwest flank appear to be levees, like those that form beside subaerial lava flows. Between them and farther downslope are two bathymetric lobes, with the uppermost one backstepping in offlap, as if formed by a later smaller flow. The summit is surrounded by crater rims. The central crater has lost little material to reworking through wave base. Many red bombs rest on the upper slopes (Figs. 7a, 7b, 7c, and 7e), much as they fell, despite the Holocene transgression. The upper section of proximal cores (Figs. 1d and 8) is a bioturbated lag built on the last tephra.

The cone’s north flank has a smooth gentle slope with a north northwest arcuate ridge presenting as a radial fissure or eroded dyke. Its central prominence shows disaggregated lava blocks and flanking areas with oxidized bombs and scoria of decreasing size, as per camera transect (Table S2A, Fig. 7a) and grabs 142 and 144 (Table S2B), typical of subaerial fissure vents.

ROV video tracked northwest from the hydraulically fluted moraine and the southwest flank of the cone, traversing its breached amphitheatre, levees, and flows. Gravel to cobble-sized scoria clasts and larger bombs and blocks are encrusted by benthos. The southwest flank has a head scarp descending to a levee-bounded, lobate pahoehoe flow. These lobes are likely hawaiite or mugearite, based on morphology (Figs. 3 and 5) and tephra compositions in core 93B. An eroded ledge (Fig. 6d) marks the outer limit of the upper flow, downlapping onto an earlier one, which reached farther south. Flow lobes extend downslope to a depth of 250 mbsl where onlapping glaciomarine deposits (Fig. 6a) contain contemporaneous IRD, delicate glassy airfall tephra, and benthic foraminifera for ¹⁴C dating (Table 1, Core 124 Fig. 8). Core 93B recovered glaciomarine sediment (light grey, stiff, dense somewhat layered silts and diamicton, and covered by ∼10 cm of soft olive-coloured sandy mud with abundant scoria, sandy tephra, shell fragments), and sand sized benthic foraminiferal tests (Fig. 8). Delicate and diverse airfall tephra are dispersed throughout the fine-grained glaciomarine mud rather than forming distinct layers. We interpret this to represent very small eruptions persisting for 2000 years. This onlapping relationship demonstrates that volcanic activity spanned local deglaciation (2.9 ka duration, Table 1) (when the seabed was subaerially exposed down to 200 mbsl atop the rebounded moraine) through later transgression and floating ice shelf conditions, when dropstones and glacial rock flour deposited from a calving ice front farther up Finlayson Channel.

In this section, we present new CHIRP seismic profiles and interpret how they inform the stratigraphy and geomorphology from before the beginning of volcanic activity on MC through its drowning and extinction. The CHIRP seismic profile through MC (Fig. 6b) shows a symmetric steep cone down to 200 mbsl. It rests on the previously formed and eroded Finlayson Moraine (Fig. 6a). Irregular bathymetry may be due to some slumping near the southwestern base of the cone, particularly near waypoint 2 just off the southwest face of the cone and waypoint 5 at about 250 mbsl. There are several knickpoints and benches along the profile of the cone. On the southwest, the volcano is mantled by glaciomarine sediment up to 131 mbsl. The gentle slopes on the terraces between 250 and 150 mbsl are about 5° (Fig. 5) like beaches. These resemble wave cut platforms and may mark stillstands during sea level rise.

Glaciomarine sediments mantle and onlap from 150–90 mbsl on the northeast face down to lower than 200 mbsl. An amphitheatre on the southwest summit, bounded by levees on both sides (Figs. 1d and 3), likely marks a breach in the cone and flows of lava or pyroclastic materials built up from about 100–150 mbsl. The CHIRP profile (Fig. 6b) shows no acoustic penetration of the cinder cone (waypoints 4–5) or eroded moraine (waypoints >5). Hard acoustic returns at positions between waypoints 2–4 correspond to gravelly lags atop the onlapping, and now eroded, glaciomarine diamictons. Southwest of position 2, the glaciomarine reflectors dip to the southwest and are folded and somewhat disrupted, as is the underlying moraine top reflector. Across the drowned continental shelf, glaciomarine reflectors are horizontal (due to deposition via ice-distal bottom-hugging turbidity currents), but here, they dip downslope and may be disrupted by deglacial motion on the Principe Laredo Fault. The Holocene hemipelagic muds show onlap towards the northeast onto the glaciomarine and earlier moraine materials (southwest of position 1 Fig. 6a and 2 Fig. 6b). The glaciomarine reflectors have apparently been disrupted (offset) by shallow neoglacial faulting. Both the glaciomarine turbidites and the moraine have hard acoustic returns indicating a lag seabed with abundant gravel armour. North northeast of position 4 (Fig. 6a), Finlayson Channel has been hydraulically overdeepened and eroded down to more than 650 mbsl. There is minor accumulation of hemipelagic mud below 500 m, but poor acoustic penetration due to a high gas content.

CHIRP seismic profile (Fig. 6c) obliquely crosses Higgins Moraine. The line break in the profile is the turn against the rocky shore of Swindle Island. A series of North northeast gently dipping, downstepping terraces skirt Price Island, are underlain by glaciomarine sediments, atop the bedrock-cored moraine (Fig. 9). Northwest along Higgins Passage, postglacial hemipelagic mud is confined to the deepest basin greater than 180 mbsl. On the profile, dashed red lines denote possible faults.

Off Higgins Moraine, in Milbanke Sound, the seabed is very flat and deeper than 250 mbsl. Elsewhere along the west coast, this deep mud accumulation marks the final facies formed during marine transgression (rising systems tract) and maximum flooding stage with near modern sea levels by 6–7 ka. The mud has a faint internal reflector and downlap or offlap midway through its nearly transparent section. Regionally this is recognized as the Holocene discordance. By 6 ka, sea levels within 2 m of present established modern current patterns and tidal residuals, so mud drifts ceased migrating and built up uniformly with steady coastline positions. Postglacial mud comprises sandy silts to silty clays with low acoustic backscatter, acoustic transparency, continuous, and draping with weak or moderate internally stratified reflections in seismic profiles. We interpret the flatness of the seabed is due to a former wave base or otherwise tractive surface, and only a few metres of postglacial Holocene hemipelagic mud, as dated in core 123 (Fig. 8).

The KH samples are located in Fig. 1c and Table S1 (“17-BJB” series of Bednarski and Hamilton 2019). The MC samples (Fig. 1d, Table S2B) were taken by clamshell grabs from the cone where the seabed retained abundant volcanic scoria, blocks, and bombs. Petrographic observations (Table S3, Figs. S2 and S3, and Supplementary File 4) delineate the petrogenesis of these volcanoes as dry, within-plate, and low volume partial melts (about 1%) from asthenosphere in the garnet stability field. Volcanic activity and deglaciation were contemporaneous from glass, IRD, and epifauna in the same samples. The eruptive products of both volcanoes are diverse, even within single facies from multiple small eruptions, rather than from a single monogenetic eruption. The petrography is inconsistent with MSV being a recurrence of the Miocene rift setting or a volcanic arc.

Both volcanoes have olivine and plagioclase with three different textures. The freshest olivine basalt and hawaiite scoria bombs on MC resemble those at KH in their glassy vesicular porphyritic texture, olivine preservation, and crystallization sequence: olivine–plagioclase–augite–magnetite. Xenocryst olivine and plagioclase are polycrystalline, fractured, eroded, embayed, resorbed, and rimmed along with rare xenocrysts of cumulate textured, exsolved, abraded accidental orthopyroxene, implying residence in a mantle lithospheric or deep crustal magma chamber, and contamination rather than direct ascent from the asthenosphere. Euhedral olivine, plagioclase, and augite phenocrysts crystallized during ascent. Both KH and MC contain undercooled glasses with hopper crystals of olivine and acicular laths of plagioclase with terminal prongs, reflecting quench in glacial waters or frigid air, especially in KH’s basal subglacial tuya facies. Irregular complex vesicles at MC indicate slower degassing with some compaction when it was subaerial.

The tephra is co-deposited with IRD and glacial rock flour silt, showing the final volcanic activity occurred during glacial retreat with floating ice in the fjords when rising sea level drowned MC. MC’s delicate glassy tephra (Fig. 10), less than 3.5 mm across, is locally derived due to lack of abrasion or alteration. Decreasing size away from the vent, and stretched or flattened vesicles, demonstrate an airfall origin. Like whole rock samples, MC tephra manifest other textures: trachytic basalts (mugearites) with flow aligned plagioclase microlites and sparsely porphyritic basalts with larger olivine, plagioclase, and rare clinopyroxene or magnetite phenocrysts. MC shows a diversity of products and additional evolution compared to KH.

MSV whole rock samples and tephra are enriched Ocean Island Basalt (OIB) sodic alkalic basalts, hawaiites, and mugearites. The data are presented in Fig. 11, Tables S4 and S5, and Supplementary File 4. MC rocks plot in the silica-undersaturated alkaline fields along the divide between alkali olivine basalt and hawaiite (Cox et al. 1979). KH rocks are all hawaiites, save one mugearite. For this level of silica, Al₂O₃ values are below the range expected for arc magmas or Hi–Al basalts. This precludes any arc setting and confirms their within-plate, alkaline character, and evolution via plagioclase and olivine fractionation rather than substantial crustal feldspar contamination (either at depth or by glacial materials).

MSV samples on the alkali-iron–magnesium (AFM) diagram (Fig. 11b) are more iron rich than the calcalkaline–tholeiite dividing line (Irvine and Baragar 1971) or than precursory Masset T-MORBs (Transitional Mid-Ocean Ridge Basalt) precluding arc or rift settings. The most primitive compositions for the MSV are more evolved, with higher iron (lower Mg numbers) and alkali and trace element enrichment, than the latest Miocene lavas to come up in this same region (Hamilton and Dostal 2001, 1993; Woodsworth 1991); thus, they are more alkaline and derived by a lower amount of partial melting. Marginally greater iron enrichment for KH and Helmet peak than MC indicates heterogeneous source, and small volume partial melts. The MSV have Mg numbers too low to be directly derived from the mantle and require some residence time and early deep fractionation of Mg-rich olivine and Ca–Al rich plagioclase prior to eruption.

MC glassy tephra show two compositions (Figs. 11 and S4; Table S5). Both cores 93B and 124 contain whole rocks for the more primitive olivine basalt through hawaiite scoria, pillows, and flows from the MSV (Table S4). Hawaiites follow a simple line of descent from alkali basalts (Figs. 11a and 11b) via low pressure fractionation of about 4:1 plagioclase:olivine, without much crustal contamination as per their petrography (Table S3, Fig. S4). Some tephra grains in core 93B are more differentiated mugearites (55% SiO₂), indicating a different ascent history and likely an additional shallower crustal magma chamber.

The trace element (TE) spidergram (Fig. 11c) and rare earth plots (Fig. 11d) for MSV are flat when normalized to the OIB values (Sun and McDonough 1989). They demonstrate MSV’s enriched OIB and within-plate character and rule out arc or Hi–Al basalt or T-MORB settings. The HREE increase (Fig. 11d) requires some residual garnet, implying either a deep source for the slight decompression melting or a shallower garnet stability field. These decompression melts are simple, within-plate, small volume, and dry, but they are locally heterogeneous both within and between MSV centres. They are mantle-derived but subsequently fractionated and differentiated at depth (as per their low MgO contents) prior to eruption. They belong to sodic–alkalic series and are as enriched as the most enriched OIBs globally. They show only a little crustal contamination by intermediate to felsic plutonic rocks, which are their main substrate here in the Alexander terrane. For additional discussion of TEs, comparisons with other OIBs, mantle sources and melting, see Supplementary File 4.

No suitable material for high precision ¹⁴C dating has been found from the limited exposures of the volcano-glacial facies of MSV above sea level and these volcanoes are too young for conventional dating techniques (^39/40Ar, fission track) to be useful. Fortunately, sea level was rising during MC’s eruption, and datable foraminifera were deposited contemporaneously with airfall tephra in the surrounding sediments. We even observed foraminifera that grew in volcanic vesicles (Fig. 12, image 27).

The benthic foraminifers from Milbanke Sound (Fig. 12) were identified using the nomenclature registered in the World Foraminiferal Database (Hayward et al. 2023), and the assemblages allow us to constrain the depth at which the sediments were deposited in the cores 093B, 123, and 124 (Figs. 1d and 8). These are based on Holocene distributional patterns on the shelf and slope of the Gulf of Alaska, which provide bathymetric ranges of faunal depth zonation (Echols and Armentrout 1980), supplemented by depths of surveys reported by Todd and Low (1967) and Bergen and O’Neil (1979).

The faunal assemblages are listed in Table S6. Core 093B (38–40 cm) is dominated by Inner Sublittoral Zone (<50 m) species. Station 124 samples (134–136 cm and 14–16 cm) contain an Outer Sublittoral Zone (155–210 m) fauna. A few of the Inner Sublittoral Zone species were also recovered here, but we are basing our depth determination on the deepest-dwelling species (Bandy and Arnal 1960; Jorissen et al. 1994; McGann 2014) as downslope transport of material is possible at the toe of the moraine or across the drowned slopes of MC. Finally, core 123 samples (217–219 and 8–10 cm) are the deepest encountered in this study (>210 m). The assemblage at 8–10 cm appears to represent a slightly deeper environment than that at 217–219 cm.

Collectively, this pattern of environments and their stratigraphic positions (and relative ages) suggests a transgressive sequence where the older core 124 has a shallower assemblage at 134–136 cm than at 14–16 cm, showing transgression and rising sea levels during glaciomarine deposition. The shallowest assemblage (∼50 m) occurs in the shallowest (131 mbsl) and most proximal core position, 93B. The deepest assemblages are in the youngest core 123, after MC finished erupting, agreeing with modern core positions for maximum flooding stage and high stand conditions in late Holocene for Milbanke Sound.

The median ages and 95% confidence limits of the calibrated ¹⁴C likelihood functions are reported in Table 1.

The RSL for Milbanke Sound (Fig. 13a) shows local crustal position with respect to sea level and time. The right side (Fig. 13b) shows the relative placement of volcanic facies and geomorphology along with elevations or depths that constrain the RSL model. The 95% confidence intervals of the calibrated ¹⁴C ages are plotted along the horizontal axis of Fig. 13a, at the vertical positions of their present depths. For the five foraminiferal assemblages described in sufficient detail to make depth determinations, the possible ranges for the paleo-sea-surface are plotted as vertical grey bars. The water depths of the cores and the inferred depths of the foraminiferal assemblages along with their ages provide constraints on the RSL. Notably, the sea surface rose from at least 14 cal ka BP, but it was well below the far-field global eustatic sea level (Lambeck et al. 2014) during the MC eruptions. Equivalently, local land levels were higher from rapid rebound and incomplete forebulge collapse, while MC continued erupting tephra as its lower portions were being draped by glaciomarine sediment.

As ice withdrew, Helmet Peak erupted beneath the thinning ice sheet to form its tuya volcano from near modern sea level to 285 masl (Fig. 1e). We interpret Souther’s (1990) description of metre size blocks of granodiorite in hyaloclastite there to indicate volcano-glacial facies. Furthermore, we interpret Helmet Peak’s northeast–southwest, ice modified form (Bednarski and Hamilton 2019), and the erosional separation from its outliers along both sides of Matheison Channel (Fig. 3) to have occurred at the time of local deglaciation and locally depressed crust prior to much glacial isostatic rebound (GIR). RSL was around 75 masl for some time before complete ice retreat to erode the Milbanke Strandflat smoothing regions of subglacial incision, while building the Swindle Island strand line and fluvioglacial outwash delta. Their incision and erosion happened during rebound-driven regression.

Darvill et al. (2022) measured ¹⁰Be exposure ages on erratics and bedrock on Aristazabal Island (northwest 50 km along strike of the Principe Laredo Fault) and Hunter Island (50 km southeast). They determined that ice retreat from this coastal line occurred at 16.0–15.1 and 16.5–15.6 ka, respectively (interquartile range). Our independent analysis infers that the ice sheet withdrew from Milbanke Sound to the Coast Mountains circa 16 ka.

With ice withdrawal, local rebound was sufficient and rapid enough to uplift the floor of Milbanke Sound, currently at ∼265 mbsl, to near sea level prior to 14 ka. This required about 300 m of uplift to build terminal moraines across Higgins passage and Finlayson Channel, as ice retreated from Milbanke Sound. Deglaciation’s last phase held floating fjord ice that grounded on the recessional moraines like the one across the mouth of Finlayson Channel. Hydraulic fluting across Finlayson Moraine (Figs. 2b, 3, 5, and 9) and overdeepening was likely caused by outburst floods from ice and water farther up Finlayson Channel.

As rebound finished, MC built up subaerially atop Finlayson Moraine. MC erupted co-mingling its ash with transgressive distal glaciomarine deposits over the span of 2.9 ± 1.6 ka (2σ) from 14.1 to 11.2 cal ka BP whereon rebound and volcanic activity ceased. Finally, the deep seafloor of Milbanke Sound became blanketed with hemipelagic mud, while shallower areas of stronger currents were winnowed and grew sponge reefs.

Following the analysis of James et al. (2000, 2005), we model the crustal response as the sum of two exponential terms. Unlike the RSL model obtained at Effingham Inlet on western Vancouver Island (Dallimore et al. 2008), the directions of the vertical motion, for the two exponential decays, at Milbanke Sound oppose each other. The first has rapid rebound on local faults from the subjacent low-viscosity upper mantle to accommodate ice unloading in the fashion of an elevator. The second, slower, later response originates from deeper mantle lateral flow related to deglaciation of the Coast Mountains (and the rest of the Cordillera) collapsing the forebulge. Under Milbanke Sound, the initial rebound is vertical asthenospheric flow from the local ice load removal, while the later forebulge collapse is lateral asthenospheric flow away from this region back towards the east to infill beneath the distally rising unloaded Cordillera.

We model the RSL with deglaciation occurring as a geodynamic impulse at t₀ = 16 ka. No model fits the constraints with younger t₀; however, t₀ = 16.5 ka works as well. Regional constraints for deglaciation make older values for t₀ incommensurate. Along with local crustal motions and regional asthenospheric flow, eustatic sea level rose through the latest Wisconsinan and continued through the Holocene. The best fit models for Milbanke Sound converge with the eustatic curve by about 7 ka.

The RSL model helps explain some puzzling geomorphic differences between KH and MC. Because KH formed earlier, it experienced profound shoreface and wave base erosion in more extensive shallow seas prior to much rebound. MC, by contrast, formed after rebound and was rapidly drowned due to forebulge collapse and eustatic sea level rise. This rapid drowning in Milbanke Sound’s deep-sheltered environment preserved MC’s form.

The CIS load addition and removal performed spatially and temporally to induce crustal motions, mantle flow, and decompression melting, which produced small volume deglacial volcanoes. We present the sequence of events across Milbanke Sound (Fig. 14) and a model showing how this applied across the Cordillera at the LGM (Fig. 15). The MSV erupted during the late Wisconsinan as the CIS withdrew from the continental shelf and retreated into the Coast Mountains. We interpret the relationship between deglaciation and volcanism around Milbanke Sound as causal rather than merely coincidental. There is some evidence globally for glacial–volcanic activity related to changing loads and stresses in the crust (Nakada and Yokose 1992; Riehle et al. 1992; Watt et al. 2013; Conway et al. 2015, 2023; Maccaferri et al. 2015; Praetorius et al. 2016; Mora and Tassara 2019).

Wilson and Russell (2020) presented a potential mechanism for deglacial volcanism and pumping caused by the CIS loading and unloading, given a background level of melting and volcanic activity in the Garibaldi-Cascades Arc. There, the degree of glacial loading and unloading and the strain rates exceed background tectonic stresses by more than two orders of magnitude (Thorson 2000; Brillon et al. 2018) and profoundly influenced shallow crustal motions and neoglacial faulting as seen in Puget Sound (Thorson 1996). The ice loads were comparable or greater farther north across the Coast Mountains.

In Milbanke Sound, there is no background melt production or volcanic activity, nor has there been since the Miocene cessation of rifting in QCB. Milbanke Sound resides on top of the eastern inverted (uplifted) edge of QCB, with its thinned crust (∼20 km;  Rohr and Dietrich 1992; Yuan et al. 1992; Hyndman and Hamilton 1993; Lowe and Dehler 1995; Flück et al. 2003), high heat flow (Lewis et al. 1991), and shallow warm asthenosphere (Hyndman 2017) left over from the rifting through 17 Ma (Yorath and Hyndman 1983).

During the Wisconsinan CIS advance, the increased load over the Coast Mountains depressed the crust, expelling the asthenosphere to flow laterally and upwards to rise up beneath Milbanke Sound and QCS. This asthenospheric flow created the forebulge and instituted decompression melting beneath that region. The glacially induced uplift rates in the range of 1 m/year resemble hotspot OIB volcanoes as modelled by Elkins et al. (2008). The difference is that here the motion is not convectively driven by a mantle plume but rather induced by changing ice loads. Due to the high heat flow and tectonic history here, the melt zone in the garnet stability field resembles that of ocean islands (Wilbold and Stracke 2006) and produced similar OIBs.

At the LGM through 18 cal ka BP (Fig. 14), the Insular Belt experienced maximal ice extent, thickness, and loading (Dyke et al. 2003; Dyke 2004; Bednarski and Hamilton 2019). Our observations in Milbanke Sound include volcano-glacial facies built on glaciated bedrock or recessional moraines, and an RSL curve based on lithofacies on Swindle Island and direct dating of foraminifers and the depth constraints of their assemblages in glaciomarine sediments with airfall tephra from the drowned flanks of MC. From these data, we reconstruct the ice retreat and sea level history as a series of cross sections from QCS to the edge of the Coast Mountains. By 18 cal ka BP, the CIS over Milbanke Sound was over 1200 m thick, creating a local glacial isostatic depression of about 300 m that exceeded the far-field global eustatic lowering of sea levels (∼120 m). Swindle Island and Milbanke Sound were covered by the CIS and depressed to below sea level at LGM.

By 16–15 cal ka BP, the Wisconsinan was ending, global eustatic sea levels were rising, and the western edge of the CIS began to float and break up. The crust below QCS began to rise due to GIR. Additionally, the thick CIS farther inland loaded the crust and generated or maintained a forebulge (Josenhans et al. 1993, 1997) across to the outer continental margin. GIR in Milbanke Sound uplifted the crust, reactivated the Principe Laredo Fault, and formed the subglacial tuya beds at KH. As the CIS thinned and withdrew back into the Coast Mountains, KH emerged from the ice to become an island volcano, while rebound and forebulge exceeded sea level rise.

The late Wisconsinan retreat (Fig. 14, 16–10 cal ka BP) and loss of more than a kilometre of ice caused crustal unloading. Strain rates greatly exceeded background tectonic values, returning a tensional state to the upper crust. We hypothesize that this final uplift and decompression event caused more melting, resupplied magma chambers, and triggered eruptions of the MSV. Magmas upwelled along the weak and ready pathway of the neoglacially reactivated Principe Laredo Fault. The ice load removal changed stress to tensional at the top of the crust (convex up) and along older faults, opening them to initiate volcanism. Simultaneously the lower crust was put into compression (concave at the bottom), and this might have additionally exerted forces to squeeze and empty the lower or subcrustal magma chambers, similar to the “pumping” actions as conceived by Wilson and Russell (2020).

By 14–10 cal ka BP, the CIS withdrew farther NE into the Coast Mountains (Clague 1985; Letham et al. 2016), leaving only local alpine glaciers in the Outer Fjord Ranges. The removal of the local ice load not only caused GIR, but potentially caused decompression melting in the mantle asthenosphere in an analogous manner to the present-day interplay of GIR and enhanced melting beneath Iceland’s Vatnajökull Ice Cap (Schmidt et al. 2013). These rising melts were able to recharge, and to reheat the subcrustal magma chambers accounting for the resorption of primary phenocrysts and xenocrysts, resetting the olivine and plagioclase compositions as observed petrographically.

After local GIR ceased and as the forebulge here collapsed, the local RSL began to rise (Fig. 14, 14–10 cal ka BP) as recorded by the progressive onlap of glaciomarine sediments and the upwards deepening of foraminiferal facies. The shallowest portion of the MC pyroclastic cinder cone (above 130 mbsl) continued to grow and cast tephra into the surrounding seas as its base was drowned. The surrounding sea likely supplied water as additional propellant for the explosive pyroclastic volcanism, until activity ceased and the MSV became dormant. Following this cessation of local crustal motions, the far-field eustatic rise took over, drowning Milbanke Sound and MC and once more returned sea level to the foot of KH.

The existence of and slower response of the forebulge is essential to understanding how the ice load and its removal led to small isolated volcanic eruptions as the ice age ended. The forebulge is observed from northern Vancouver Island to Baranof Island, a region of about 800 × 150 km² (Shugar et al. 2014). In contrast to the inner limit of the forebulge of Shugar et al. (2014), in our interpretation, Milbanke Sound is profoundly affected by the forebulge, though it lies only 40 km from Shearwater (near Bella Bella), which shows no ¹⁴C evidence of synglacial uplift at 13.4 cal ka BP from two dates on marine material currently terrestrial. Thus, we consider the forebulge, at least locally, to be fault-bound, and to have behaved like an elevator rather than a long lever-arm tilting from a hinge line. The forebulge collapse completed a few hundred years earlier on Goose Bank ∼100 km to the south southwest (Luternauer et al. 1989; Barrie et al. 1991). The final stage (Fig. 14, 10–0 cal ka BP) shows the modern DEM profile and sea level, after the complete loss of the CIS and cessation of ice-induced crustal motions.

The maximal ice load at LGM is shown as a Cordilleran cross section (Fig. 15). The ice load formed on the highest part of the Coast Mountains depressing the crust. This put the upper crust into compression and expelled some asthenosphere to flow laterally and upwards beneath the thinned crust of the Insular Belt and former edge of QCB. A few kilometers of relief on the lithosphere asthenosphere boundary (LAB) caused a small amount of decompression melting (1%) in the garnet stability field below, during upwards flow of up to 1 m/year as per the model for global OIB’s (Elkins et al. 2008). This mantle upflow allowed both a small amount of melting and melt separation. By LGM, the base of the lithosphere under Milbanke Sound was in tension and magmas rose to fill chambers, where differentiation occurred via fractional crystallization of olivine and plagioclase. Ice load removal returned the local crustal stress profile to its normal interglacial condition, releasing somewhat differentiated, but still mantle-derived OIB-type sodic alkalic basalts and hawaiites. Positive feedback occurred here because glacial isostatic unloading also allowed the lithosphere to locally rise, producing more decompression melting, rejuvenating the deep magma chambers, and generating the complex phenocryst xenocryst assemblages seen petrographically. Ultimately there was only a small amount of decompression, so only small local mantle heterogeneities partially melted and small volume volcanic products were delivered over a span of about 5 ka to the MSV. The Principe Laredo Fault and ancillary faults acted as relief valve conduits for deglacial volcanism. This model explains the small volumes of volcanic materials and multiple eruptions across the span of a few thousand years as the last ice age ended.

While all volcanoes have their own controlling conditions involving sources, magma chambers, conduits, tectonic settings, local faults, and other details, we suggest that perhaps some of the neoglacial volcanic activity across the Canadian Cordillera (Grove 1974; Hickson and Edwards 2001) shares similarities to the interpretation here in MSV. Most of the Canadian Cordillera has been glaciated and deglaciated in the Wisconsinan (Dyke 2004). Most of the Cordillera has a relatively shallow and warm LAB (Hyndman and Canil 2021). Unlike the Garibaldi-Cascades Arc, there is no evidence of background magma flux beneath much of the Cordillera, yet there are abundant small cinder cones and volcano-glacial facies along its length. Given the very low tectonic strain rates and motions well inside the continental margin, the application and removal of more than a kilometer of ice load within a few thousand years is likely the largest and most changeable set of conditions the crust has seen. In the model (Fig. 15), we suggest other neoglacial volcanic products may have a similar origin to the MSV via ice-induced asthenospheric upflow, decompression melting, and volcanic activity concentrated at the time of ice load removal. If the ice load similarly reactivated crustally penetrative faults such as the Queen Charlotte Transform, renewed activity on Bowie Seamount (Cousens et al. 1985; Cousens 2010) during the Wisconsinan eustatic lowstand might have a similar cause. Ice load application and removal on faults along the edges or within the Intermontane Belt might similarly have caused small deglacial cinder cones along its edges such as Castle Rock, Atlin Volcanic Field, Alligator Lake in the Northern Cordilleran Volcanic Province (Edwards and Russell 1999), Wells Gray (Hickson and Vigouroux 2014), and Summit Lake (Francis et al. 2010). Undoubtedly there are other significant factors such as mantle hot spots, local heterogeneities, and longer within-plate records of local volcanic activity to explain widespread but small volume late Quaternary volcanoes across the Cordillera. Nonetheless, the immense working of thin lithosphere and a warm asthenosphere by a thick ice sheet is a powerful process that should be considered for magma genesis, particularly for the small volume and widespread volcano-glacial facies of the Canadian Cordillera.

Milbanke Sound, along BC’s central coast, features a 43 km line of deglacial volcanoes, distributed along the southern limit of the Principe Laredo Fault, a former strike–slip boundary within QCB. The fault was reactivated as a conduit for deglacial magmatism and a zone of weakness to accommodate stresses arising from glacial loading and unloading. The erosional morphologies in Milbanke Sound, both above and below sea level, together with the volcanic and sedimentary facies provide the basis for understanding the geological history since the end of the last ice age and provide the constraints to understand how these small volcanoes formed.

The volcanoes were subglacial as the ice thinned, or deglacial and subaerial as the ice cap withdrew. This setting includes the originally subaerial MC, which currently lies 200–43 mbsl. Airfall tephra deposited with glaciomarine sediments around the volcano. Benthic foraminifera in these sediments date from 14.1 to 11.2 cal ka BP, and their depth ranges provide a unique setting to determine an RSL curve.

The volcanoes are petrologically and geochemically OIB, with the implication that they result from very slight decompression melting less than a few percent and subsequent differentiation within a deep magma chamber. Crustal rebound following deglaciation was sufficient to change crustal stresses and to release the deep magmas, providing the causal mechanism for volcanism to follow deglaciation in this setting.

We modelled the local RSL as a superposition of the rapid drop due to postglacial rebound starting at 16 ka, plus a slower longer range forebulge collapse starting at 15 ka eventually meeting the rising global eustatic sea level. Our “ice then fire” model may be applicable across the Canadian Cordillera, suggesting that crustal load of ice sheets and their removal can account for some of the recent small volume igneous activity even for within-plate settings.

KH lies in the unceded territory of the Kitasoo/Xai’Xais First Nation (Kitasoo Spirit Bear Conservancy), and the field work was done under the auspices of the Kitasoo Band Office. The authors wish to thank Vernon Brown for guiding and student assistants Jacintha Brown, Robbie Duncan, Tina Lobbes, and Mercedes Robinson-Neasloss. LIDAR survey was partly supported by the Canada Foundation for Innovation and NSERC grants to Brian Menounos, University of Northern British Columbia. Thanks are due to Robert Kung of the GSC for his assembly of the DEM.

The authors would like to thank the captain and crew of Canadian Coast Guard Ship Vector for supporting the sampling programs in Milbanke Sound. We thank Tammy Norgard and the Department of Fisheries and Oceans for the use of ROV data on McGregor Bank. We would also like to thank Annika Meijer and Felix Parkinson for help with preliminary seabed mapping work.

Radiocarbon measurements were made at the National Ocean Sciences Accelerator Mass Spectrometry facility of the Woods Hole Oceanographic Institute in Massachusetts, USA.

We are particularly thankful for the stimulating discussions and encouragement from our late colleague Roy Hyndman. We are grateful for helpful reviews from Rene Barendregt and Chris Conway.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the Canadian or US Governments.

All primary research data are included in the manuscript and Supplementary material.

Conceptualization: TSH, RJE, JMB, CDS

Data curation: TSH, RJE, ZL, JMB, CDS, MLM, BJLJ

Formal analysis: TSH, RJE, ZL, JMB, CDS, MLM, BJLJ

Funding acquisition: JMB, CDS

Investigation: TSH, RJE, ZL, JMB, CDS, MLM, BJLJ

Methodology: TSH, RJE, JMB, CDS, MLM, BJLJ

Project administration: RJE, JMB, CDS

Resources: CDS

Software: RJE, ZL

Supervision: TSH, RJE, JMB, CDS

Validation: TSH, RJE, JMB

Visualization: TSH, RJE, ZL, JMB, CDS, MLM, BJLJ

Writing – original draft: TSH, RJE

Writing – review & editing: TSH, RJE, ZL, BJLJ

Funding for this research was provided by NRCan’s Marine Geoscience for Marine Spatial Planning. Multibeam bathymetry was acquired by the Canadian Hydrographic Service and provided to the Geological Survey of Canada through CHS MOU No. 2022-0324-1260-N. The foraminiferal and radiocarbon portions of this study were funded by the Natural Hazards, Coastal and Marine Hazards and Resources Program of the U.S. Geological Survey Pacific Coastal and Marine Science Center. This is NRCan contribution number 20230229.

Supplementary data are available with the article at https://doi.org/10.1139/cjes-2023-0080.