Dissecting 20 million years of deep-water forearc sediment routing using an integrated basin-wide Bayesian chronostratigraphic framework

Forearc basins are complex and diverse geodynamic systems in which the subduction of oceanic crust exerts a strong control on sediment supply (Dickinson and Seely, 1979; Ridgway et al., 2011; Hessler and Fildani, 2015), subsidence (Moxon and Graham, 1987; von Huene and Scholl, 1991; Fildani et al., 2008; Stern, 2011), magmatism (e.g., Moxon and Graham, 1987; Kortyna et al., 2014; Cecil et al., 2018), and the structural evolution of the margin (Dickinson and Seely, 1979; Cashman et al., 1992; Noda, 2016; Bishop et al., 2019). Submarine channel systems in forearc regions facilitate the transport of sediment from terrestrial sources to oceanic trenches (Klaus and Taylor, 1991; Bourget et al., 2011; Fig. 1). Their sedimentary deposits can record the signal of geodynamic processes in both the terrestrial sediment source area (e.g., underplating, uplift, and magmatism) and the basin (e.g., trench accretion and subduction erosion). However, isolating these signals in ancient forearc strata is difficult, in part, because forearc basin fills are commonly tectonically deformed or fragmented both during and after basin sedimentation (Mitchell et al., 2010; Takano et al., 2013; Hessler and Fildani, 2015; Orme et al., 2015; Noda, 2016; Greene and Surpless, 2017; Herriott et al., 2019). The submarine slopes of forearc regions are often characterized by structural complexity, which can lead to local controls on sediment thickness, timing of sedimentation, and sediment routing system morphology (Bourget et al., 2011; Paquet et al., 2011; Takano et al., 2013; Hessler and Fildani, 2015; Orme and Graham, 2018; Fig. 1). As such, careful correlations and widespread analyses may be needed to disentangle forearc margin evolution from deep-water sedimentary deposits.

Resolving the chronology of sedimentary basin fills, including submarine channel deposits, is paramount for understanding sediment routing system evolution and establishing links to tectonic processes (Kimbrough et al., 2001; Fildani et al., 2009; Takano et al., 2013; Daniels et al., 2018; Englert et al., 2018; Pinter et al., 2018; Herriott et al., 2019). This requires the creation of chronostratigraphic frameworks that accurately constrain the architecture of the deposits (e.g., deposit geometry, spatial distribution, and paleoflow direction) as well as the timing of sedimentation along large outcrop belts (Fildani et al., 2008; Takano et al., 2013). When successful, the evolution of deep-water sediment delivery systems can be used to elucidate the rates and spatial extents of subtle geodynamic processes (Daniels et al., 2018; Johnstone et al., 2019; Englert et al., 2020). Basin frameworks have traditionally been constructed through field mapping and the compilation of biostratigraphic, magnetostratigraphic, tephrachronologic, and U-Pb detrital zircon geochronologic datasets (Kimbrough et al., 2001; Fildani et al., 2009; Eddy et al., 2017; Pinter et al., 2018; Huang et al., 2022). These methods have limitations that cannot be addressed in a robust fashion when used individually as the sole age constraints across a study area. Statistical sampling methods for chronostratigraphic data (Johnstone et al., 2019; Halverson et al., 2022) provide an avenue for integrating different chronometers and decreasing the impact of method-specific limitations; however, so far they have only been applied in single 1-D stratigraphic sections through intervals of interest. The extension of these methods to 2-D basin-wide cross sections or faulted 3-D exposures where the distribution of chronometers (biostratigraphic, magnetostratigraphic, and detrital zircon ages) may vary could significantly improve methods for capturing and measuring large, basin-scale processes.

To better understand the long-term spatio-temporal controls on sedimentation in forearc settings, we studied the protracted record of deep-water sediment transfer in outcrops of the upper Nanaimo Group, Nanaimo Basin, Canada. These outcrops form a 160-km-long depositional strike-oriented belt composed of numerous submarine channel-system deposits. We combined 2199 m of measured stratigraphic section, mapping of major sediment routing-system deposits, and an integrated chronologic dataset composed of detrital zircon maximum depositional ages, biostratigraphic zones, and magnetostratigraphic ages from the upper Nanaimo Group to: (1) for the first time, create a detailed chronostratigraphic framework of a basin margin using a Bayesian numerical model; (2) use this integrated approach to assess the accuracy of detrital zircon maximum depositional ages; and (3) use it to interpret the controls on submarine channel-system evolution in the Nanaimo Basin. Our novel chronostratigraphic approach provides unprecedented insight into along-strike, deep-water sediment routing and forearc basin evolution while demonstrating the application of a statistical technique to a basin-scale chronostratigraphic dataset. We also include considerations that should be accounted for when assembling large chronostratigraphic models.

Subduction of oceanic crust along the western margin of North America from the late Paleozoic through Cenozoic resulted in the sequential docking of terranes, the generation of magmatic arcs, and ultimately, the formation of the North American Cordillera (see Monger et al., 1982; Umhoefer and Miller, 1996; Cecil et al., 2018). Subduction of the Farallon Plate was roughly normal to the Pacific margin of North America during the Mesozoic until the early Late Cretaceous, when the Kula Plate broke away from the Farallon Plate along the Kula-Farallon ridge. This resulted in the Kula Plate trending northward from the Kula-Farallon ridge, and the Farallon plate tending southward from the ridge (Woods and Davies, 1982). Spreading at the Kula-Farallon ridge occurred as early as 84 Ma (Woods and Davies, 1982; Engebretson et al., 1985) and was continuous by 79 Ma (Seton et al., 2012). Rifting of the Farallon plate rotated the subduction vector of the Kula Plate to be oriented highly oblique to the margin (i.e., subducting to the north-northeast; Woods and Davies, 1982; Doubrovine and Tarduno, 2008).

Magmatism associated with this subduction is recorded by the Coast Mountains Batholith, which is the eroded remnant of a Mesozoic continental magmatic arc on the western margin of Canada (Monger et al., 1982; Rusmore and Woodsworth, 1991; Miller et al., 2009; Cecil et al., 2018, 2021; Fig. 2A). Magmatism in the Coast Mountains Batholith spanned the Jurassic through Eocene, with the locus migrating eastward at a rate of ~2.5 km/m.y. starting at 120 Ma due to thickening of the continental arc (Friedman et al., 1995; Gehrels et al., 2009; Cecil et al., 2018). As a result, no magmatic activity occurred on the arc's western margin, nearest to the Nanaimo Basin, throughout much of the basin's history (Santonian to Danian).

The Nanaimo Basin formed atop the Wrangellia terrane during the Late Cretaceous after Wrangellia had docked to western North America and was sutured by the Coast Mountains Batholith (Monger et al., 1982; Mustard, 1994; Monger and Price, 2000; Tikoff et al., 2022; Fig. 2A). The basin is situated in a forearc position, between the Coast Mountains Batholith, at the basin's eastern limit, and the Kula/Farallon subduction zone, which is located on the western edge of Wrangellia (Muller and Jeletzky, 1970). The southern margin of the Nanaimo Basin is in thrust contact with the northwest Cascade thrust system (Fig. 2A), which is composed of imbricated allochthonous terranes of diverse age and composition (Brandon et al., 1988; Brown, 2012). Terranes of the northwest Cascade thrust system were thrust in place at ca. 84 Ma, and sediment was shed into a restricted portion of the southernmost Nanaimo Basin that is directly adjacent to the thrust system (Brandon et al., 1988; Coutts et al., 2020).

Evolution of the Nanaimo Basin was closely tied to the subduction of the Kula/Farallon Plate beneath the North America Plate (Englert et al., 2018; Kent et al., 2020; Coutts et al., 2020). Changes in subduction obliquity are correlated with enhanced subsidence of the Nanaimo Basin, the initiation of basin-wide deep-water deposition, and changes in sediment source areas connected to the Nanaimo Basin (Englert et al., 2018; Coutts et al., 2020). Much of the basin's fill was deeply eroded (England and Calon, 1991; England et al., 1997), and the subduction complex was rifted from the margin by the Kula Plate during the Eocene (Garver and Davidson, 2015; Fig. 2A). Nanaimo Group formations presently crop out as part of the Cowichan fold-and-thrust belt, which formed at ca. 50 Ma (England and Calon, 1991). There are minor structures in the study area, and dips predominantly range from 30° to 5° and generally dip to the east (Fig. 2).

The Nanaimo Basin is filled by the Nanaimo Group, which is divided into informal lower and upper portions based on regional lithostratigraphic correlations and interpreted depositional settings (Muller and Jeletzky, 1970; Mustard, 1994; Figs. 2A and 2B). The lower Nanaimo Group consists of both terrestrial and marine strata that crop out on the eastern side of Vancouver Island and in limited outcrops on the southern Gulf Islands of western British Columbia (Mustard, 1994; Jones et al., 2018; Kent et al., 2020; Figs. 2A and 2B). The distribution of lower Nanaimo Group deposits was influenced by topography on the Wrangellian basement, notably the Nanoose and Cowichan basement uplifts (Muller and Jeletzky, 1970; Pacht, 1984; Mustard, 1994; Kent et al., 2020; Fig. 2A). The upper Nanaimo Group, the focus of this study, consists of deep-water deposits, including three coarse-grained formations (DeCourcy, Geoffrey, and Gabriola formations) that comprise the majority of exposed bedrock on the Gulf Islands as well as minor outcrops on Vancouver Island (Figs. 2A and 2B). These units primarily consist of large-scale, channel-form–shaped, sandstone and conglomerate deposits (1400–4100 m wide and 60–280 m thick) that are interpreted as long-lived slope-channel systems that crosscut the margin of the Nanaimo Basin between the Late Cretaceous and Paleocene (England and Hiscott, 1992; Rowe et al., 2002; Bain and Hubbard, 2016; Englert et al., 2018, 2020). The intervening fine-grained formations (Cedar District, Northumberland, and Spray formations) occur lateral to and between coarse-grained lithostratigraphic units and are interpreted to represent slope strata, including proximal and distal overbank deposits (Bain and Hubbard, 2016; Englert et al., 2018; Figs. 2B–2E). The deposition of coarse-grained formations was, in part, a response to magmatism, underplating of subducted rocks, and subduction erosion, which affected the basin and hinterland (Katnick and Mustard, 2003; Englert et al., 2018; Coutts et al., 2020).

The age of the Nanaimo Group was constrained by biostratigraphy (McGugan, 1962, 1979, 1982; Muller and Jeletzky, 1970; Ward, 1978; Haggart and Ward, 1989; Galloway et al., 2023; Ward et al., 2012; McLachlan et al., 2018, 2019), paleomagnetic analysis (Enkin et al., 2001; Ward et al., 2012), and detrital zircon maximum depositional ages (MDAs; Englert et al., 2018, 2020; Huang et al., 2019, 2022; Kent et al., 2020; Coutts et al., 2020). Micro- and macro-fossils retrieved from the fine-grained slope strata were used to subdivide the Nanaimo Group into biostratigraphic zones that range from Turonian to Maastrichtian (McGugan, 1982; Haggart and Ward, 1989; McLachlan et al., 2018, 2019); however, biostratigraphic ages cannot be determined for non-fossiliferous, coarse-grained formations (e.g., Gabriola Formation on Hornby Island; Muller and Jeletzky, 1970; Ward et al., 2012). Due to the close proximity of the Nanaimo Basin to the Coast Mountains Batholith, MDAs derived from abundant near-depositional age detrital zircon provide accurate constraints on the depositional ages of coarse-grained strata (e.g., Englert et al., 2018; Kent et al., 2020). Recent analysis by Englert et al. (2020) and Huang et al. (2022) demonstrated the compatibility of ages derived from biostratigraphic, magnetostratigraphic, and detrital zircon analyses.

The provenance of the Nanaimo Group is constrained by petrologic point-counting datasets (e.g., Pacht, 1984; Mustard, 1994; Coutts et al., 2020) and U-Pb detrital zircon (e.g., Matthews et al., 2017; Mahoney et al., 2021). Basin-wide petrological results generally show a change from transitional arc-derived to basement uplift-derived sediment (sensu Dickinson et al., 1983) throughout the deposition of the lower and upper Nanaimo Group. Samples proximal to the northwest Cascade thrust system exhibit chert-rich and lithic frameworks (e.g., transitional recycled) sourced from the obducted chert units and exotic terranes in the thrusts (Pacht, 1984; Mustard, 1994; Coutts et al., 2020). U-Pb detrital zircon data indicate there are four dominant zircon facies within the Nanaimo Group. Zircon facies 1 is derived from crystalline rocks with major age modes at 150 Ma and 90 Ma, which is consistent with documented periods of high-flux magmatism in the Coast Mountains Batholith (Gehrels et al., 2009; Miller et al., 2009; Cecil et al., 2018; Matthews et al., 2017; Huang et al., 2019; Coutts et al., 2020) and influence from the Wrangellian basement. Minor sediment sources for the lower Nanaimo Group included Wrangellian basement (scattered age modes of ca. 450 Ma and 350 Ma). Zircon facies 1 is found in the lower Nanaimo Group and older formations of the upper Nanaimo Group. Detrital zircon facies 2 is found in the upper portions of the upper Nanaimo Group and is derived predominantly from the Coast Mountains Batholith, with major modes between 80 Ma and 70 Ma and minor age modes from 250 Ma to 110 Ma. Unlike much of detrital zircon facies 1, sediment in the detrital zircon facies 2 was also derived from a source area that contained Proterozoic zircon with major age modes near 1700 Ma and 1380 Ma (Matthews et al., 2017; Coutts et al., 2020; Isava et al., 2021; Mahoney et al., 2021). The source of Proterozoic zircon in the Nanaimo Group is contentious, with either the Mojave Sonoran region of California or the Lemhi Subbasin of Idaho being proposed (Matthews et al., 2017; Mahoney et al., 2021). Depth-profiling and cathodoluminescence imaging of Nanaimo Group zircon reveals that some Cretaceous and all Proterozoic detrital zircon exhibit textures that are indicative of Late Cretaceous (100–72 Ma) high-grade metamorphism of the sediment source area (Boivin et al., 2022). These results were interpreted to indicate that a significant proportion of detrital zircon in the upper Nanaimo Group was derived from a Cretaceous-aged, high-grade metamorphic rock that was rapidly uplifted between 84 Ma and 72 Ma and recycled into the basin (Boivin et al., 2022). Detrital zircon facies 3, 4, and 5 are restricted to the southernmost portion of the basin, proximal to the northwest Cascade thrust system in the lower Nanaimo Group. These facies received sediment sourced from various terranes of the northwest Cascade thrust system and the Coast Mountains Batholith. Detrital zircon facies 3 and 4 have zircon derived from the Turtleback and Yellow Aster terranes (discontinuous populations from 2550 Ma to 1100 Ma and from 600 Ma to 300 Ma; Brown, 2012; Coutts et al., 2020) along with some Coast Mountains Batholith. Detrital zircon facies 5 is nearly entirely sourced from the Easton Metamorphic Complex (single age mode of 150 Ma).

Outcropping strata of the upper Nanaimo Group were examined on Gabriola, Mudge, Valdes, Galiano, Mayne, Saturna, and Vancouver islands, and observations were integrated with data collected by Bain and Hubbard (2016; Hornby and Denman islands) and Englert et al. (2018; Gabriola, Galiano, and Mayne islands; Fig. 2). One hundred stratigraphic sections were measured, totaling 2199 m of centimeter-scale descriptions. Grain size, bed thickness, sedimentary structures, and bedding contact geometry were used to define facies and interpret depositional processes (e.g., Miall, 1996). Facies associations represent a succession of rock (meters to tens of meters in thickness and tens to hundreds of meters in length/width) defined by a range of genetically related facies. The outcropping facies associations were mapped using a differential GPS (dGPS) with submeter resolution. Stratigraphic components are mapped depositional bodies (tens to hundreds of meters thick and hundreds to thousands of meters wide) that are defined by their internal facies associations, overall geometry (e.g., width, thickness, and basal surface character), and lateral facies relationships. Lastly, channel-system deposits are composed of multiple stratigraphic components that represent the deposits of long-lived sediment routing systems. The orientations of stratigraphic components were interpreted from measurements of paleoflow direction (n = 4207) that were made from erosional scours, sole marks (flute, tool, and groove casts), imbricated clasts, and cross-stratification.

Field data, as well as the top bedrock lithology present in publicly available water well records (Government of British Columbia, 2015), were projected onto a two-dimensional plane oriented perpendicular to bedding (method of Englert et al., 2018). This allowed mapped stratigraphy to be viewed perpendicular to the dip of the beds (i.e., where structural tilt was removed). Because paleoflow was generally westward while bedding dips to the east, strata were largely viewed in a depositional strike-oriented cross section using this method. The stratigraphic architecture of upper Nanaimo Group deposits was interpreted on these projections to demonstrate relationships between coarse- and fine-grained deep-water deposits.

Appropriate uncertainty propagation during detrital zircon data generation is critical to the determination of MDAs due to the reliance on individual grain uncertainties in MDA calculation algorithms. Over or underestimation of uncertainty can lead to calculated MDAs that are too old or too young, respectively, particularly in datasets with continuous distributions of zircon dates down to the depositional age of the sample. Importantly, session-to-session variation in the magnitude of excess variance (ε), a large component of random uncertainty that is calculated through statistical analysis of the dispersion of a relatively small number (eight or nine) of measurements of a validation reference material measured with the unknowns (Horstwood et al., 2016; Matthews and Guest, 2017), can result in fictitious variations in MDA.

In this study, we combined detrital zircon U-Pb data first reported by Matthews et al. (2017), Englert et al. (2018), and Coutts et al. (2020). These data were acquired over ~48 measurement sessions using the method of Matthews and Guest (2017). To address the issue of session-to-session variations in ε, and produce the most consistent and reproducible MDAs possible, we modified the uncertainty propagation approach of Matthews and Guest (2017) to harness the large number of measurements of validation reference materials in the combined dataset.

Our uncertainty propagation followed a multi-sessional approach. Excess variance and average Pb signal intensity was calculated for each of three to five validation references in each measurement session. These data were combined from all ~48 sessions to produce a single calibration curve between ε and Pb signal intensity for each of the isotopic ratios (²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U). The Pb signal intensity of the unknown was then used to calculate an appropriate ε value for each datum.

The use of a single calibration curve across all 48 sessions assures that uncertainty propagation is not creating spurious variations in calculated MDAs. MDAs from upper Nanaimo Group strata exposed on Hornby and Denman islands were calculated using data propagated in a similar manner and agree with biostratigraphic ages that constrain the fine-grained units of the basin (e.g., Englert et al., 2020).

Maximum depositional ages were calculated using the youngest grain cluster at 2σ uncertainty (YGC 2σ) method of Dickinson and Gehrels (2009). The YGC 2σ method was demonstrated by Coutts et al. (2019) to rarely calculate MDAs younger than the true depositional age (<<1% of cases). When combined with high-throughput acquisition methods (e.g., Matthews and Guest, 2017; Pullen et al., 2014), as was done in the datasets used here, the YGC 2σ method calculates MDAs with small uncertainties by including numerous grains within a weighted average. This allows for superior resolution of geologic processes with shorter durations (e.g., Daniels et al., 2018; Coutts et al., 2019). Given the long duration of the basin fill history, the YGC 2σ method also enables the basin fill to be subdivided using MDAs with small uncertainties (0.97–2.12 m.y. 2σ). Information regarding the age, uncertainty, sample location, and original publication of the data used to calculate the maximum depositional ages presented here is provided in File S1 in the Supplemental Material¹.

Maximum depositional ages calculated from detrital zircon dates can only provide an upper limit on the depositional age of sampled strata (Kimbrough et al., 2001; Nelson, 2001; Rainbird et al., 2001; Dickinson and Gehrels, 2009). Because of this, other chronostratigraphic constraints (e.g., biostratigraphic and magnetostratigraphic ages) must be employed to assess the accuracy of calculated MDAs and hone their interpretation (Johnstone et al., 2019; Halverson et al., 2022). We incorporated the MDAs compiled by Coutts et al. (2020) with local magnetostratigraphic (Ward et al., 1997; Enkin et al., 2001; Englert et al., 2020) and biostratigraphic constraints (McGugan, 1979; Ward et al., 2012) into a Bayesian Monte Carlo age model (see Johnstone et al., 2019). Following the methods of Johnstone et al. (2019), synthetic ages are drawn from age probability distributions that represent the potential age of the input samples based on the available constraints. Zircon MDAs are represented by maximum/survival functions that go to zero probability at the maximum 2σ uncertainty limit of the MDA. Biostratigraphic and magnetostratigraphic constraints are represented by a uniform probability distribution between the lower and upper age limit of each biozone or chron. In its most simplistic form, each Bayesian model is composed of stratigraphically ordered probability distributions representing a stratigraphic section or transect. In each iteration of the model, a single age is drawn from each age probability distribution. If the sequence of ages generated violates the law of superposition, the ages are rejected. If they obey the law of superposition, they are accepted and added to the depositional age distribution (t_d) for each sample. This process is repeated numerous times (2000 times in our study) to fill out the t_d. The mode of the t_d represents the best-fit distribution of each sample, given the surrounding constraints. For a more detailed explanation, see Johnstone et al. (2019).

The majority of biostratigraphic and magnetostratigraphic constraints on the age of the upper Nanaimo Group is found on Hornby, Denman, and Gabriola islands (McGugan, 1979; Ward et al., 1997; Enkin et al., 2001; Ward et al., 2012; Englert et al., 2020), with only sparse samples or broad age constraints elsewhere in the basin (McGugan, 1979; Muller and Jeletzky, 1970). Many fine-grained successions with known biostratigraphic or magnetostratigraphic constraints are lateral to coarse-grained deposits of interest here, with no directly mappable correlations (Bain and Hubbard, 2016; Englert et al., 2020). As such, biostratigraphic and magnetostratigraphic constraints were only selected if the sample's stratigraphic position could be interpreted with a high degree of certainty. Two Bayesian age models were created, one integrating data from Hornby and Denman islands with the other focused on Gabriola Island (Supplemental File S2).

To validate the MDAs against surrounding age constraints, the MDAs were compared to their associated depositional age distributions (t_d) calculated by the Bayesian model. Synthetic ages were randomly generated from Gaussian distributions fitting the parameters of the MDAs (age and 2σ uncertainty) and subtracted from the corresponding t_d distribution. This created an uncertainty distribution, which represents the difference between the t_d (the most accurate depositional age given the surrounding age constraints) and the MDA. If MDAs are accurate constraints on depositional age, this uncertainty distribution will yield a mode near zero. If the intercalated MDAs are older than the surrounding age constraints (e.g., represent the reworking of previously deposited material), there will be a large difference between MDA and t_d, and the uncertainty distribution will have a large mode. If the MDAs are too young for the surrounding age constraints (e.g., affected by lead loss or sample contamination), the model will fail to run because this scenario breaks the law of superposition. By comparing the MDA to the corresponding t_d distributions, we can assess the accuracy of the MDAs in areas where additional chronostratigraphic constraints exist and extend our confidence in the MDAs to guide interpretations in other areas of the basin where only MDAs are available and additional constraints are sparse.

Eight lithofacies were identified in the upper Nanaimo Group (L1–L8; Table 1). All facies are interpreted to have been deposited through processes associated with sediment gravity flows. Thick-bedded, amalgamated conglomerate and sandstone (L1); thick-bedded, amalgamated sandstone (L2); thick-bedded, non-amalgamated sandstone (L3); and discontinuous, thinly interbedded sandstone and siltstone (L4) are interpreted to be deposited through traction and suspension processes associated with energetic, low- to high-density turbidity currents (i.e., Ta, Tb, and Tc divisions of Bouma, 1962; R and S divisions of Lowe, 1982; T_b1–T_b4 associations of Postma and Cartigny, 2014; Table 1). Thin-bedded, fine-grained sandstone and mudstone (L5) is attributed to traction and suspension processes of low-density turbidity currents (i.e., Tb through Te divisions of Bouma, 1962; Table 1). Massively bedded mudstone and siltstone (L6) is interpreted to be deposited by dilute turbidity currents and through suspension settling of mud and silt-sized sediment (i.e., Te divisions of Bouma, 1962; Table 1). Mudstone-matrix conglomerates (L7) are interpreted to be deposited by cohesive debris flows, while disorganized blocks of sandstone and siltstone (L8) are interpreted to be deposited through mass-wasting events, including large slumping events (e.g., Middleton and Hampton, 1976; Lowe, 1976; Alves, 2015; Qin et al., 2017; Table 1). Six lithofacies associations (FA1–FA6) were identified based on the spatial relationships and proportions of encompassing lithofacies.

Facies association 1 (FA1) is predominantly composed of clast-supported and sandstone matrix-supported conglomerate (L1) interbedded with thick-bedded, amalgamated sandstone (L2; 60%–85%), while non-amalgamated sandstone (L3 and L5), discontinuous thin-bedded sandstone and siltstone (L4), and mass-transport deposits (L7 and L8) compose minor proportions (15%–25%, 0%–5%, and 0%–5%, respectively). FA1 deposits make up successions 5–15 m thick that are bound by undulatory to concave-up basal surfaces with <5 m of relief (Fig. 3A). At their base, FA1 successions are commonly characterized by mass-transport deposits (L7 and L8); discontinuous interbedded sandstone and siltstone (L4); thinly interbedded, fine-grained sandstone and mudstone (L5); or amalgamated conglomerate (L1). These are often overlain by a fining upward succession (L1 → L3/L5; Fig. 4A). The coarsest facies (L1 and L2) transition laterally to finer-grained facies (L2 and L3) over tens of meters.

Facies association 1 deposits are interpreted to represent the product of gravel-dominated submarine channels (e.g., Hickson and Lowe, 2002; Kane et al., 2007; Greene and Surpless, 2017). Conglomeratic deposits (L1) and amalgamated sandstones (L2) are attributed to deposition in channel thalwegs, where repeated high-energy flows scoured the substrate (Cronin et al., 2000; Kane et al., 2007; Hubbard et al., 2008; Kneller et al., 2020). Amalgamated to non-amalgamated sandstone beds (L3 and L5) reflect deposition adjacent to the thalweg, where flows were less energetic, which led to increased preservation of entire flow units (Mutti and Normark, 1987; Hickson and Lowe, 2002; Jobe et al., 2010; Hubbard et al., 2014).

Thick-bedded sandstone (L2 and L3) comprises the majority of facies association 2 (FA2; 50%–60%), while thinly interbedded sandstone and siltstone (L5); discontinuous thinly interbedded sandstone and siltstone (L4); mass-transport deposits (L7 and L8); and thick-bedded, amalgamated sandstone and conglomerate (L1) compose smaller proportions (15%–25%, 0%–5%, 0%–5%, and <1% respectively). FA2 deposits form successions 5–10 m thick that have undulatory to concave-up basal surfaces with <5 m of observed relief (Fig. 3B). At their bases, FA2 successions are characterized by amalgamated sandstone, thin successions of thinly interbedded sandstone and siltstone (L5), or mass-transport deposits (L7 and L8), which are overlain by a generally fining upward succession from L2 to L3/L5 (Fig. 4B). Amalgamated deposits (L2) commonly transition laterally into non-amalgamated and then finer-grained facies (L3 and L5) over tens of meters. In some localities, nested wedges of sandstone (L2) 5–20 m long and 75–200 cm thick dominate FA2 successions (Fig. 3C). These wedges stack opposite to the dominant paleoflow direction (i.e., migrate upstream).

Facies association 2 deposits are interpreted to be the fills of sand-dominated submarine channels due to the abundance of erosional surfaces filled with high-density turbidites (e.g., Mutti and Normark, 1987; Macauley and Hubbard, 2013; Casciano et al., 2019). Amalgamated sandstones (L2) are attributed to channel thalweg processes, whereas non-amalgamated deposits (L3 and L5) reflect processes that are active away from the thalweg, where finer-grained sediments are more likely to be deposited and preserved (e.g., Walker, 1966; Macauley and Hubbard, 2013; Jobe et al., 2017; Bell et al., 2018). Nested wedges of sandstone were recently interpreted as the reworked deposits of supercritical flow bedforms that are observed within or in association with submarine channels (Hage et al., 2018; Englert et al., 2021).

Facies association 3 (FA3) is mainly composed of thick-bedded sandstone (L2 and L3) (90%–100%), while thin-bedded, fine-grained sandstone and siltstone (L5) compose minor proportions (0%–10%). FA3 deposits form successions 4–10 m thick that are bound by basal surfaces that are concordant to underlying strata (Fig. 3D). The majority of FA3 successions are characterized by uniform bed thickness, as well as similar maximum grain size and bed amalgamation. However, examples of bed fining and thinning or bed coarsening and thickening are locally observed (Figs. 3D and 4C). FA3 successions are laterally continuous at the scale of outcrop exposures (50–900 m) and show minimal lateral facies transitions.

The abundance of coarse-grained turbidites (L2 and L3) that exhibit de-watering and traction structures suggest rapid deposition from high-density turbidity currents (Lowe, 1982; Mutti and Normark, 1987; Talling et al., 2012; Lowe et al., 2019; Fig. 4C); however, a lack of regular lateral facies transitions and the tabular nature of the beds are indicative of non-channelized deposition (e.g., Normark, 1970; Shultz and Hubbard, 2005; Spychala et al., 2015). Bed thickness trends and changes in sandstone bed amalgamation are related to the evolution of weakly confined flows on intraslope lobes (e.g., Pyles, 2008; Prélat and Hodgson, 2013; Spychala et al., 2015).

Facies association 4 (FA4) is composed primarily of L4 (70%–85%), while thick-bedded sandstone (L2 and L3), massively bedded siltstone (L6), mass-transport deposits (L8), and discontinuous thinly interbedded sandstone and siltstone (L4) compose smaller proportions (0%–5%, 15%–30%, 0%–5%, and 0%–2% respectively; Fig. 3F). FA4 deposits form sedimentary successions 5–20 m thick with basal surfaces that are concordant with the underlying strata. FA4 successions thin and fine-upward from L3 to L6 or have uniform bed thickness and grain-size distributions (Fig. 4D). L4 deposits form discrete intervals within FA4 successions that are 0.75–2 m thick and composed of numerous discontinuous sandstone beds. L4 sandstone beds range from flat lying to inclined (dipping toward or opposite to the regional paleoflow direction) and are truncated by surrounding mudstone. FA4 is laterally continuous at the scale of outcrop exposures (5–50 m; Fig. 3E), except where it is juxtaposed laterally with deposits of FA1 and FA2.

Facies association 4 is interpreted to represent proximal out-of-channel deposition, in environments that include levees (e.g., Normark et al., 2002; Kane and Hodgson, 2011; Morris et al., 2014; Hansen et al., 2015), splays (e.g., Wild et al., 2005; Cross et al., 2009; Lowe et al., 2019), and terraces (Hansen et al., 2015, 2017). The high abundance of low-density turbidites (L5) capped with massive mudstone and siltstone (L6) is consistent with low-energy, turbulent flows in an out-of-channel environment. Rare thick-bedded, coarse-grained sandstone (L2 and L3) with dewatering structures is indicative of rapid overbank splay deposition (Hubbard et al., 2009; Lowe et al., 2019). Discontinuous thinly interbedded sandstone and siltstone packages (L4) are interpreted as the record of sediment bypass and bedforms formed under supercritical flow conditions (e.g., Postma and Cartigny, 2014; Bain and Hubbard, 2016; Cornard and Pickering, 2019; Stacey et al., 2019), which are attributed to flows escaping channel confinement (e.g., Fildani et al., 2006). Mass-wasting deposits (L8) result from slumps that mobilized blocks of fine-grained material along levee slopes both toward and away from the related channel (e.g., Deptuck et al., 2003; Hubbard et al., 2009).

Mass-transport deposit facies (L7 and L8) compose >95% of facies association 5 (FA5), while other facies (L1–L6) compose a minor proportion (<5%; Fig. 3G). FA5 deposits compose 5- to 10-m-thick sedimentary successions (Fig. 4E), with concordant to undulatory bases that exhibit <5 m relief; however, basal surfaces of FA5 deposits are often poorly expressed in outcrop, and the deposits could be thicker. FA5 successions do not exhibit any variation vertically or laterally at the scale of outcrop exposures (5–50 m wide; Fig. 3G).

Facies association 5 (FA5) is attributed to submarine mass-wasting of material from farther upslope of the strata that were examined. FA5 deposits are often related to erosion and undercutting during the initial inception of submarine conduits (Deptuck et al., 2003; Mayall et al., 2006), slumping of levee material into and away from the channel (e.g., Deptuck et al., 2003), or mass failures on the continental slope (e.g., Saller and Dharmasamadhi, 2012; Alves, 2015).

Massively bedded siltstone (L6) comprises the majority of facies association 6 (FA6; 75%–85%), while thick-bedded sandstone (L2 and L3), interbedded sandstone and siltstone (L5), and mass-transport deposits (L8) compose a small proportion (5%–10%, 15%–25%, and 0%–5%, respectively). FA6 comprises successions 5 m to >50 m thick and has basal surfaces that are concordant with the underlying facies associations. FA6 successions have few vertical trends, with sandstone beds (L2 and L3) distributed randomly throughout the succession (Figs. 3H and 4F). FA6 is laterally continuous at the scale of outcrop exposures (5–50 m wide).

FA6 is interpreted to represent relatively quiescent deposition in distal levee, splay, or background slope depositional settings (Mutti, 1977; Wild et al., 2005; Di Celma et al., 2011; Hansen et al., 2017). Thick-bedded, massive to laminated mudstones are interpreted to be deposited through suspension settling of silt and mud on the continental slope (Bouma, 1962; Hill, 1984; Kane et al., 2007; Newport et al., 2018; Boulesteix et al., 2019). Thin-bedded sandstones are interpreted as the deposits of unconfined, low-density turbidity currents consistent with distal out-of-channel deposition (e.g., Mutti and Normark, 1987; Kane et al., 2007; Hansen et al., 2017). Rare, thick-bedded, coarse-grained sandstones are interpreted as the deposits of high-density turbidity currents, which we attribute to distal splay processes (Wild et al., 2005; Cross et al., 2009; Lowe et al., 2019).

The facies associations of the upper Nanaimo Group generally comprise three distinct channel-form–shaped stratigraphic components (CSCs), as well as three tabular stratigraphic components (TSCs).

CSC1s have channel-form geometry and range from 50 m to 170 m thick and are 2–4 km wide. They are composed of both conglomerate-dominated and sandstone-dominated submarine channel deposits (FA1 and FA2), as well as minor thin-bedded deposits (FA4) and mass-transport deposits (FA5; Figs. 5A and 6A). Although the internal architecture of CSC1s is difficult to discern due to high degrees of amalgamation, they locally consist of smaller scale (15- to 30-m-thick and >100-m-wide) sedimentary bodies that are also channel form in shape (Fig. 7). Examples of CSC1s are exposed along Dinner Point (DeCourcy Formation, Mayne Island), East Point (Geoffrey Formation, Saturna Island), and on Mount Galiano (Geoffrey Formation, Galiano Island; Fig. 2E).

CSC1s are interpreted to be the deposits of long-lived, gravel-dominated submarine conduits (e.g., Morris et al., 1989; Hickson and Lowe, 2002; Kane et al., 2007; Hubbard et al., 2008; Kneller et al., 2020). The presence of internal channel-form bodies suggests that CSC1 deposits represent composite submarine channel units, which are equivalent to channel complexes or channel complex sets (Sprauge et al., 2002; Beaubouef, 2004; Di Celma et al., 2011; Macauley and Hubbard, 2013; Cullis et al., 2018).

Low-aspect-ratio, sandstone-dominated channel units are channel form in geometry and range from 50 m to 140 m thick and are 2–4 km wide. CSC2s are primarily composed of sandstone-dominated channel deposits (FA2) and minor thin-bedded (FA4) and mass-wasting (FA5) deposits (Figs. 5B and 6A). CSC2s are composed of sedimentary bodies that are 15–30 m thick and ~150 m wide, which are commonly difficult to discern due to significant amalgamation (Figs. 3B and 6A). Examples of CSC2s are exposed at Monarch Head (DeCourcy Formation, Saturna Island), Duke Point (DeCourcy Formation, Vancouver Island), and East Point (Geoffrey Formation, Saturna Island; Figs. 2D and 2E).

CSC2s are interpreted to result from processes of erosion, sediment bypass, and deposition in long-lived submarine conduits (e.g., Mutti and Normark, 1987; McHargue et al., 2011; Englert et al., 2020). The presence of internal channel-form bodies suggests that CSC2 deposits are highly composite submarine-channel units that are equivalent to channel complexes or channel-complex sets (Sprauge et al., 2002; Di Celma et al., 2011; Macauley and Hubbard, 2013; Cullis et al., 2018; Casciano et al., 2019).

CSC3s are 30–120 m thick, 3–7 km wide, and composed of sandstone-dominated channel deposits (FA2) with minor mass-transport deposits (FA5; Figs. 5C and 6). The internal architecture of CSC3s is difficult to observe in many instances due to bed amalgamation (Figs. 5C and 7). The highly amalgamated deposits are generally characterized by little lateral change in facies associations. Where exposed, internal channel-form bodies are ~5–15 m thick and >400 m wide (Figs. 5C and 7). Examples of these architectural components are exposed at Taylor Point (DeCourcy Formation, Saturna Island), East Point (Geoffrey Formation, Saturna Island), and Dragon Bay (Gabriola Formation, Gabriola Island; Fig. 2E). The internal architecture of CSC3s can also be composed of FA2 dominated by L2 beds that are discontinuous and lenticular. Examples of this expression of CSC3s are exposed at Taylor Point (DeCourcy Formation, Saturna Island), Twin Beaches Peninsula (Gabriola Formation, Gabriola Island), and Edith Point (Gabriola Formation, Mayne Island; Figs. 2D and 2E).

CSC3s are interpreted to record sediment transfer within broad channel belts (e.g., Clark and Pickering, 1996; Kane et al., 2007; Hubbard et al., 2008; Bayliss and Pickering, 2015; Pinter et al., 2018; Butler et al., 2020). Internally, thin channel forms are interpreted as the deposits of weakly confined channels (e.g., Moody et al., 2012; Brunt et al., 2013; Bayliss and Pickering, 2015; Pemberton et al., 2016). Laterally offset stacking of these elements is consistent with laterally migrating channels prone to avulsions within a broad fairway (e.g., Clark and Pickering, 1996; Moody et al., 2012; Bayliss and Pickering, 2015). Thick successions of FA2 composed of discontinuous beds of L2 have been attributed to segments of broad conduits characterized by fields of supercritical flow bedforms (Cornard and Pickering, 2019; Englert et al., 2021).

Sandstone-dominated tabular components are 30–100 m thick, at least 0.2–1 km wide, and composed primarily of FA3 (Fig. 5D). FA3s stack to form thick and laterally extensive TSC1 successions. Paleoflow measurements of TSC1s have similar mean paleoflow azimuths (i.e., NNW). Examples of TSC1s are exposed along Vineyard Road (Gabriola Formation, Galiano Island) and the southwestern cliffs of Valdes Island (Gabriola Formation; Figs. 2D and 2E).

Amalgamated FA3s of TSC1s are interpreted to form thick successions of weakly confined depositional lobes (e.g., Mutti and Normark, 1987; Prélat and Hodgson, 2013; Pinter et al., 2018; Lowe et al., 2019; Shultz and Hubbard, 2005). Although deposit tabularity is reflective of largely unconfined depositional processes, consistent unidirectional paleocurrent measurements within individual stratigraphic components and the vertical stacking of successive FA3 deposits with little heterogeneity is consistent with the presence of some degree of lateral topographic confinement (e.g., Takano et al., 2005; Pinter et al., 2018).

TSC2s are 25–250 m thick and at least 100 m wide, with measurements limited by the extent of intertidal outcrop (Fig. 6A). They are primarily composed of thinly interbedded siltstone and sandstone (FA4), less frequent siltstone successions (FA6; Fig. 5E), and rare mass-transport deposits (FA5; Fig. 3F). FA4 deposits within TSC2s commonly show upward vertical bed-thinning trends in 5- to 20-m-thick packages (Figs. 3E and 3F). Mass transport deposits (FA5) are dominantly composed of rotated blocks (L8) of the surrounding TSC2 strata. TSC2s are commonly observed adjacent to coarse-grained channel-form units (CSC1–3; Fig. 6A), show evidence for syn-sedimentary normal faulting locally (Figs. 5E and 5F), and exhibit little lateral heterogeneity over outcrop extents (Fig. 6A). Examples of TSC2s are exposed at Veruna Bay (DeCourcy Formation, Saturna Island) and Miners Bay (Spray Formation, Mayne Island; Fig. 2E).

TSC2s are interpreted to be associated with out-of-channel deposition due to the tabular geometry of their deposits, along with their occurrence lateral to channel-form units (CSC1–3; e.g., Kane and Hodgson, 2011; Hansen et al., 2015). Component facies associations, internal geometry, and upward bed-thinning trends are consistent with levee successions documented in modern and ancient channel systems (Manley et al., 1997; Cronin et al., 2000; Kane et al., 2007; Normark et al., 2002; Deptuck et al., 2003; Morris et al., 2014). Within TSC2s, FA4 deposits are interpreted to represent proximal areas of levees, FA6 deposits are associated with distal areas of levees, and FA5 deposits are interpreted to represent the slumping or mass wasting of levee deposits along the flanks of channels (Kane et al., 2007).

TSC3s are 25 m to >250 m thick and 30 m to >200 m wide, composed predominantly of siltstone successions (FA6), with minor thin-bedded sandstone and siltstone (FA4), and mass-transport deposits (FA5; Figs. 5E and 5F). TSC3s exhibit little lateral heterogeneity and are predominantly observed between coarse-grained deposits (CSC1–3; Fig. 6). Examples of TSC3 deposits are exposed at Murder Point (Cedar District Formation, Saturna Island) and False Narrows (Northumberland Formation, Gabriola Island; Figs. 2D and 2E).

TSC3s are interpreted to be deposited by hemipelagic settling processes active on the slope of the Nanaimo Basin during periods of non-channelized deposition (e.g., Hill, 1984) or on the distal portions of levees that bound channel systems. Thick FA6 successions are representative of fine-grained, suspension-settling slope processes (Bouma, 1962; Newport et al., 2018; Boulesteix et al., 2019). Isolated thin-bedded, out-of-channel successions (FA4) are interpreted to represent distal splay deposits (e.g., Lowe et al., 2019), while mass-wasting (FA5) in the continental slope environment is commonly driven by slope instability (e.g., Saller and Dharmasamadhi, 2012; Alves, 2015).

Coarse-grained stratigraphic components (i.e., CSC1–3 and TSC1) stack and amalgamate to form the deposits of long-lived submarine channel systems (Fig. 8). These deposits strongly control the outcrop extent of lithostratigraphic formations of the upper Nanaimo Group due to differential erosion of the coarse-grained and fine-grained facies associations (Figs. 2C–2E and 6). Coarse-grained lithostratigraphic formations are composed of channel-system deposits, which make up topographic highs and cliffs on the interior of islands that are often highly discontinuous along depositional strike due to their channelized nature (Fig. 6). Fine-grained lithostratigraphic formations are dominantly composed of out-of-channel and background slope deposits of TSC2 and TSC3, which are exposed along island edges (Fig. 6). The presence of fine-grained deposits predominantly along the edges of islands suggests that much of the strata eroded between islands along the transect are related to fine-grained deposition (Fig. 2A; Katnick and Mustard, 2003; Bain and Hubbard, 2016).

Each coarse-grained lithostratigraphic formation of the upper Nanaimo Group is composed of numerous channel-system deposits (Fig. 8). The DeCourcy Formation is composed of at least six distinct channel systems, including both conglomerate-dominated and sandstone-dominated channel-form stratigraphic components (CSC1–3; Fig. 8). Channel systems of the Geoffrey Formation are typically composed of conglomerate-dominated and sandstone-dominated stratigraphic components (CSC1 and CSC2) flanked by thick successions of heterolithic-dominated tabular components (TSC2; Fig. 8). Channel systems of the Gabriola Formation in the northern portion of the outcrop belt are composed of conglomerate-dominated and sandstone-dominated, low-aspect-ratio channel forms (CSC1 and CSC2); however, the formation is composed of high-aspect-ratio, sandstone-dominated channel forms (CSC3) and sandstone-dominated tabular components (TCS1) in the southern portion of the outcrop belt (Fig. 8).

Since the architecture of the upper Nanaimo Group is complex and lithostratigraphic units are highly diachronous (Englert et al., 2020), the evolution of the upper Nanaimo Group is discussed primarily as a function of time and channel-system activity.

The Bayesian chronostratigraphic framework successfully fit depositional age distributions (t_d) to the input MDAs and surrounding biostratigraphic and magnetostratigraphic constraints (Fig. 9A). Uncertainty distributions of each sample (representing the difference between the t_d distributions and maximum depositional age) are narrow, with standard deviations ranging from 1.27 m.y. to 3.34 m.y. at 2σ (Figs. 9B–9D). The uncertainty distribution for all combined samples has a mode of −0.6 m.y. and ranges from −6.5 m.y. to 4.5 m.y. (Fig. 9D). A cumulative uncertainty distribution with a mode near 0 demonstrates that the MDAs are very similar to the t_d distribution. This indicates that detrital zircon MDAs provide a consistent constraint on the depositional timing of these deposits compared to the results of the Bayesian model, which include biostratigraphic and magnetostratigraphic datasets. Thus, minimal lag time is interpreted between zircon crystallization and deposition in the basin, and MDAs from areas in the basin without dense biostratigraphic or magnetostratigraphic control (i.e., Galiano, Mayne, and Saturna islands) can be treated as accurate representations of depositional ages.

Integration of field-mapping data with the results from the Bayesian numerical models and detrital zircon MDAs from previously published data (Englert et al., 2018; Coutts et al., 2020; Figs. 8A and 9; Supplemental File S1) enables the creation of a chronostratigraphic framework that constrains the along-strike variability in channel-system activity, including its duration (Fig. 10A). Channel-system deposits of the upper Nanaimo Group cluster spatially into three long-lived sediment fairways spaced 50 km and 80 km apart (Fig. 8). The northern fairway is composed of channel-system deposits exposed on Denman and Hornby islands that record ~16 m.y. of channel processes (Figs. 2C and 10; Englert et al., 2020). The northern fairway also records a major change in sediment provenance during the Maastrichtian (ca. 72 Ma), which has been attributed to the transition from local Coast Mountains Batholith drainages to expanded catchment areas that included areas east of the Coast Mountains Batholith (Matthews et al., 2017; Coutts et al., 2020; Isava et al., 2021; Mahoney et al., 2021) or the introduction of a new sediment source within the Coast Mountains Batholith at that time (Boivin et al., 2022). The central fairway is composed of channel-system deposits exposed across Vancouver, Mudge, Gabriola, and Valdes islands (Figs. 2D and 10), which record ~8 m.y. of sedimentation. Lastly, southern fairway channel-system deposits are exposed on Galiano, Mayne, and Saturna islands and record 18 m.y. of channel-system processes (Figs. 2E and 10). Both the central and southern fairways experienced the same change in provenance observed in the northern fairway during the Campanian (ca. 84 Ma; Coutts et al., 2020; Boivin et al., 2022).

Upper Nanaimo Group strata were divided into six time intervals that represent different phases of sedimentation within the Nanaimo Basin. In cases where the confidence intervals of MDAs within a single fairway overlap due to similar depositional ages or large uncertainties, the longevity of a single-channel system or the boundary between two-channel systems is challenging to interpret (e.g., samples 17 and 18; Fig. 10B).

Initial channelization in the upper Nanaimo Group is recorded in the southern fairway by two successive channel systems composed of high-aspect-ratio channel forms (CSC3) followed by a channel system composed of low-aspect-ratio, sandstone-dominated channel forms (Fig. 10B). These channel-system deposits have west-directed paleoflow directions that are perpendicular to the continental margin. However, individual channel-form deposits record paleoflow directions ranging from northwest to southeast. (Fig. 10C).

During time interval 2, the northern fairway was active and characterized by southwest-directed sediment bypass, incision, and limited deposition (Fig. 10; see Englert et al., 2020). These processes are recorded by TSC2 deposits, which are exposed on the western coastline of Denman Island (Bain and Hubbard, 2016; Englert et al., 2020; Figs. 2C, 7B, and 10). The southern fairway was characterized by the presence of a southwest-directed channel system that incised >300 m prior to filling with intercalated CSC1 and CSC2 (Figs. 8B and 10). This large-scale incision is interpreted due to the higher stratigraphic position of an older channel system (i.e., DeCourcy Formation, Saturna Island MDA number 3 in Fig. 10) as compared to that of a younger channel system (i.e., DeCourcy Formation, Mayne Island; MDA numbers 6 and 7 in Fig. 10). As there are no stratigraphic relationships exposed between islands, the MDAs and relative stratigraphic positions are used to inform this relationship.

During time interval 3, the northern fairway was occupied by a southwest-directed channel system composed of numerous CSC1s (Fig. 10). The central fairway was occupied by a northwest-directed channel system composed of CSC3, and the southern fairway records the continued presence of a southwest-directed channel system, including nested CSC3s (Figs. 8B and 10).

During time interval 4, all three fairways were occupied by gravel-dominated channel systems composed primarily of CSC1 with minor CSC2 (Englert et al., 2018; Fig. 10). Channel-system deposits range in orientation from south-directed to west-directed (Fig. 10) and are largely perpendicular to the continental margin.

During time interval 5, the northern fairway experienced a prolonged 4 m.y. phase of deep erosion (~550 m) and sediment bypass (see Englert et al., 2020; Figs. 8B and 10). The central and southern fairways were occupied by a single-channel system, which was oriented to the northwest and parallel to the continental margin. This channel system records both a significant change in sediment routing-system orientation (~90° clockwise) and in depositional style compared to the underlying channel systems (Fig. 10). This system is represented stratigraphically by a series of CSC3s that interfinger with packages of TSC1s (Fig. 8B). CSC3s are most common in the northern and southern extents of the outcrop belt (Gabriola and Mayne islands), with TSC1-dominated deposits preserved in the middle segment (Galiano and Valdes islands; Fig. 8B).

During time interval 6, the northern fairway was occupied by a westward-directed, gravel-dominated channel system composed of CSC1 and CSC2 (Fig. 10). There are no strata of this age exposed in the central or southern portions of the outcrop belt (Fig. 10D).

Stratigraphers have a great interest in constraining how sediment routing systems evolve in time and space (e.g., Hadler-Jacobsen et al., 2005; Romans et al., 2016; Sharman et al., 2018; Gradstein et al., 2020). Key deposits (e.g., basal channel deposits, large conglomeratic sequences, or the onset of low deposition rates) are often correlated to geodynamic drivers to clarify the linkages between source and sink (Kimbrough et al., 2001; Fildani et al., 2008; Englert et al., 2018). The creation of chronostratigraphic frameworks has largely relied on biostratigraphic, magnetostratigraphic, or tephrachronologic interpretations from fine-grained deposits that are intercalated between or are along-strike of coarse-grained deposits (e.g., Kimbrough et al., 2001; Pinter et al., 2018; Gradstein et al., 2020). More recently, detrital zircon MDAs allowed the stratigraphic community to create chronostratigraphic frameworks detailing the ages of coarse-grained systems (Dickinson and Gehrels, 2009; Greene and Surpless, 2017; Daniels et al., 2018; Sickmann et al., 2018). However, relying on a single age determination method makes a basin framework vulnerable to limitations of the selected method (e.g., preservation bias of ashes or fossils in fine-grained deposits). The Bayesian chronostratigraphic method developed by Johnstone et al. (2019) allows multiple types of chronometers from an area of interest to be integrated into a single framework, enabling ages to be validated and refined in many cases.

The distribution of MDAs and additional age constraints across the study area can dictate the development and application of a Bayesian model. Johnstone et al. (2019) produced a single transect (i.e., one stratigraphic section) of 13 MDAs that were optimized with two U-Pb ash ages and five fossil assemblage age constraints. The resulting t_ds were used directly for the age interpretation, as all of the additional constraints were included in the transect. In this study, MDAs were obtained across three basin transects (north, central, and southern); however, sufficient biostratigraphic and magnetostratigraphic data are only available in the northern and central transects. Bayesian models were used to calculate t_ds for MDAs in the northern and central transects but could not be generated for the southern transect. As such, the validity of the MDAs was assessed by computing the difference between the t_ds and MDAs in the basin. This approach leveraged the Bayesian model's ability to refine ages to provide best-fit ages and the large amount of additional age constraints found in the northern and central transects, while keeping the required MDAs in the southern transect.

Although numerous constraints can be incorporated into a Bayesian age model, samples must be curated to ensure that the models are accurate. When working in sedimentary sequences with complex architecture, the genetic relationships between stratigraphic units of interest with MDAs and those constrained by other methods need to be understood (e.g., crosscutting fluvial architecture, complex stacking of shoreface deposits, and incised valley fill; Fig. 11). For example, large erosional surfaces, like those with 350 m of relief that bound channel deposits in the upper Nanaimo Group (Fig. 8B), can juxtapose coarse-grained channel-fill deposits against much older deposits (Fig. 11A). In some cases, the genetic relationships between sampled units may be known from careful mapping and stratigraphic analyses; but if correlations cannot be demonstrated, then the additional constraint should not be used in the model. In the framework presented here, additional age constraints were largely incorporated from fine-grained strata between coarse-grained bodies that clearly obeyed the laws of superposition, whereas along-strike relationships between coarse-grained units of interest and biostratigraphically controlled fine-grained deposits were unknown, and these samples were removed from the analysis. Lastly, in the process of building Bayesian age models, ages may require reinterpretation or standardization. Older biostratigraphic and magnetostratigraphic ages should be updated to new biozone and chron definitions as well as the current geological time scale (Gradstein et al., 2020, at the time of this study), and compiled radiogenic ages require that all error propagations be standardized. Issues with previously published age constraints may become apparent, such as miscorrelations of magnetostratigraphic chrons due to a lack of biostratigraphic age constraints (Englert et al., 2020), MDAs that appear too old due to the reworking of older strata, or MDAs that appear unreasonably young due to sample contamination or methodological limitations of the chosen MDA method (Coutts et al., 2019). In these cases, reinterpretation of the input age constraints or their removal from the dataset may be necessary. As these considerations cover many disciplines of chronostratigraphy, it is best to assemble a multidisciplinary team to assess all interpretations and assumptions that are present in the underlying data while curating samples that will be used in the model.

The Bayesian framework used here demonstrated that MDAs are accurate chronometers throughout the upper Nanaimo Group succession. This allowed for lithostratigraphic formations (Muller and Jeletzky, 1970) to be broken up into individual component sedimentary systems dated by MDAs. Some formations of the upper Nanaimo Group represent diachronous sedimentation events, which are composed of channel systems that do not overlap in age (e.g., channel-system deposits that compose the DeCourcy Formation; Fig. 10, time intervals 1–4), while others are truly synchronous sedimentation events and composed of channel-system deposits that overlap in age (e.g., channel systems that comprise the Geoffrey Formation; Englert et al., 2018; Fig. 10, time interval 4). Resolving these intricacies was accomplished by retrieving multiple MDAs from the same lithostratigraphic unit across the entire basin rather than applying a single age (biostratigraphic or MDA) to lithostratigraphically correlated strata. Providing absolute ages for these deposits allows for regional stratigraphic models to be integrated with other datasets (e.g., hinterland thermochronology, eustatic sea-level, climate, etc.), and potentially disentanglement of the relative influences of allogenic processes on stratigraphic evolution. As there are likely many lithostratigraphic formations that exhibit similar undetected diachroneity globally, decreasing the uncertainty of age constraints in these units through widespread sampling and high-resolution dating methods will allow subtle external controls on stratigraphic architecture, regional correlations, and channel-system evolution to be recognized (Daniels et al., 2018; Coutts et al., 2019).

Forearcs undergoing oblique subduction are often compartmentalized into numerous along-strike depocenters due to the partitioning of strain into trench-normal and trench-parallel vectors (McCaffrey, 1992; Noda, 2016). These depocenters, or subbasins, are often separated by fault-bounded bathymetric highs and have different subsidence patterns (Takano et al., 2013; Berglar et al., 2017; Zhang et al., 2022). This is well documented in the Sumatra forearc basin, where highly oblique subduction segmented the forearc region into distinct subbasins separated by dextral strike-slip and transpressional structures (Berglar et al., 2017). In this Sumatra example, along-strike variations in sediment thickness between subbasins are observed (Berglar et al., 2017). Similarly, the along-strike variations in the timing of sedimentation, channel-system architecture, and large-scale channel-system reorientation identified in the upper Nanaimo Group dataset may have resulted from the differential subsidence of forearc subbasins (i.e., the fairways recognized here) driven by oblique subduction of the Kula Plate (Fig. 12). The differential subsidence of the fairways may have impacted the accommodation, slope angle, or confinement in each fairway, driving changes in sediment routing system architecture both along-strike in adjacent fairways (e.g., Fig. 10, time interval 3) or among successive channel systems in a single fairway (e.g., Fig. 7B, southern fairway). Additional processes with along-strike variability, such as drainage expansion (Matthews et al., 2017; Englert et al., 2018; Coutts et al., 2020) and hinterland deformation through underplating (Sauer et al., 2018; Coutts et al., 2020) could produce along-strike variations in sediment supply; however, the relative contributions of these processes cannot be isolated in deep time with this dataset.

Major changes in the basin, such as the initiation of deep-water sedimentation and wide-spread channel-system activity, were previously linked to changes in plate convergence, which supports the significant influence of subduction dynamics on forearc basin evolution (Englert et al., 2018; Coutts et al., 2020). The dataset presented here additionally documents a large change in sediment routing system architecture and orientation in the southern portion of the Nanaimo Basin at ca. 69 Ma (Fig. 10). This sediment routing system reorientation is not coincident with a change in sediment provenance (Coutts et al., 2020), but occurred synchronously with an incisional phase in the northern transect and an ~1 m.y. counterclockwise rotation of the subduction vector, which resulted in the Kula plate vector becoming more parallel to the margin (i.e., decreasing convergence; Doubrovine and Tarduno, 2008; Fig. 12). As such, we hypothesize that this phase of channel reorientation and incision may be related to the propagation of transpressional or transtensional structures in the basin. In modern forearc basins, this type of deformation often propagates in the direction of subduction (e.g., Malatesta et al., 2013; Berglar et al., 2017) and occurs during periods of plate reorganization (e.g., Hessler and Fildani, 2015; Malatesta et al., 2016; Andjić et al., 2018; Zhang et al., 2022).

The >60-km-long, margin-parallel channel system that persisted from ca. 68 Ma to 66 Ma (i.e., southern Gabriola Formation, Fig. 10) is characterized by a broad fairway with evidence for weak channel confinement (e.g., CSC3 and TSC1). Channel-system architecture changes along the midsection of the sediment routing system from broad channels (i.e., CSC3) to structurally confined lobe units (i.e., TSC1), which indicates a decrease in confinement or slope gradient (e.g., Takano et al., 2005; Pyles, 2008; Moody et al., 2012). Although the structure that reoriented the channel system is not observed in outcrop, the channel-system architecture is comparable to that of other structurally influenced and confined channel systems globally (Takano et al., 2005; Bayliss and Pickering, 2015; Pinter et al., 2018). These types of systems often consist of sediment routing systems with lobe-like deposits, channels, and sediment wave/bedform fields (Takano et al., 2005; Pyles, 2008; Cornard and Pickering, 2019). Structures such as active faults or mobile salt in the subsurface are widely reported to modify the slope angle and confinement along segments of deep-water sediment routing systems over distances of tens to hundreds of kilometers (e.g., Adeogba et al., 2005; Moody et al., 2012; Butler et al., 2020). Notably, these downstream changes in confinement or slope angle commonly modify channel-system architecture (Takano et al., 2005; Clark and Cartwright, 2009; Kane et al., 2010; Pinter et al., 2018). In the absence of large, mappable bounding faults, the evolution of channel systems across the multiple fairways provides indirect evidence and important context for the structural development of the westernmost Cordillera, including controls on sedimentation.

The integration of basin-scale stratigraphic analysis with a Bayesian chronostratigraphic framework based on extensive sedimentological and stratigraphic measurements, 37 detrital zircon maximum depositional ages (MDAs), biostratigraphic data, and magnetostratigraphic measurements spatially and temporally constrains the along-strike variability in submarine channel systems across a 160-km-long segment of a Late Cretaceous continental margin. The spatial characteristics (e.g., width, thickness, basin position, and paleoflow direction) and temporal characteristics (e.g., system longevity) were constrained for 12 channel-system deposits, encompassing over ~20 m.y. of sediment routing. Channel-system deposits cluster spatially within three fairways that preserve a protracted history of sediment transfer. The along-strike variations in deep-water sediment routing as well as the large-scale reorientation of channel systems in the southern portion of the basin may be related to basin segmentation and deformation of the basin due to oblique subduction of the Kula Plate. The framework presented here identifies strong links between subduction vectors and both the activity and position of sediment routing systems in forearc basins, as well as demonstrates the utility of integrated basin-wide datasets for constraining the timing of coarse-grained sediment routing system deposits in other convergent margins globally.

Funding for this research was provided by a Natural Sciences and Engineering Research Council of Canada Discovery Grant (RGPIN-2018-04223) to S. Hubbard and a Queen Elizabeth II scholarship from the University of Calgary to D. Coutts. Insightful and enjoyable field assistance was provided by Heather Gray, Ben Daniels, Sarah Southern, and Miranda Walters. The authors also thank Sam Johnstone for his help implementing his Bayesian chronostratigraphic approach. The manuscript benefited greatly from the reviews of Larissa Hansen, Devon Orme, Kelly Thomson, and an anonymous reviewer. Lastly, the authors thank the many people of the Gulf Islands who provided local knowledge, accommodation, and transportation throughout our time in the islands.