Systematic vertical organization of matrix-rich and associated matrix-poor sandstones in ancient deep-marine slope and basin-floor deposits

Deepmarine matrix-rich (>20% clay and silt content) strata that are unlike classical turbidites (sensuBouma, 1962; Lowe, 1982) or debrites have been increasingly recognized in both modern and ancient turbidite systems (e.g., Lowe and Guy, 2000; Haughton et al., 2003, 2009; Hodgson, 2009; Kane and Pontén, 2012; Talling, 2013; Terlaky and Arnott, 2014). Termed slurry beds (Lowe and Guy, 2000), linked debrites (Haughton et al., 2003), co-genetic debrite-turbidite beds (Talling et al., 2004), hybrid event beds (Haughton et al., 2009), transitional flow deposits (Kane and Pontén, 2012), and bipartite facies (Angus et al., 2019), these strata typically consist of a sandy basal part overlain sharply by a muddier, mud clast–bearing upper part, with the contact being either planar (e.g., Haughton et al., 2009; Talling, 2013; Pierce et al., 2018) or irregular (e.g., Haughton et al., 2009; Patacci et al., 2014; Fonnesu et al., 2015, 2016, 2018; Southern et al., 2015), or less commonly gradational and marked by an intervening unit of alternating light (mud-poor) and dark (mud-rich) colored bands (e.g., Barker et al., 2008; Davis et al., 2009; Pierce et al., 2018). Additionally, these composite strata are reported to occur downflow of clean (i.e., less mud-rich) sand and upflow of argillaceous (i.e., mud-rich) sand over distances of hundreds of meters to tens of kilometers in what has been interpreted to be an along-flow proximal- to- distal depositional continuum (e.g., Fonnesu et al., 2015, 2016, 2018; Kane et al., 2017; Spychala et al., 2017a; Pierce et al., 2018; Angus et al., 2019; Ningthoujam et al., 2022). The central focus of these studies has been aimed largely at understanding either the physical origin of the composite strata at the individual bed scale (e.g., Davis et al., 2009; Haughton et al., 2009; Baas et al., 2011; Talling, 2013; Fonnesu et al., 2015, 2016; Angus et al., 2019; Koo et al., 2020) or their stratigraphic and geographic distribution at the stratal element scale (several-meters-thick stratal package comprising two or more event beds and in which >50% of the beds are composed of a single sedimentary facies) (e.g., Haughton et al., 2009; Kane and Pontén, 2012; Patacci et al., 2014; Terlaky and Arnott, 2016; Kane et al., 2017; Spychala et al., 2017b; Pierce et al., 2018). Lacking, however, is an in-depth discussion of how these beds stack vertically and therefore how the sedimentological conditions that deposited them changed in time. Accordingly, this study describes the vertical stacking character and spatial pattern of matrix-rich and associated matrix-poor lithologies observed at the bed and stratal unit scales in slope and proximal basin-floor deposits of the Neoproterozoic Windermere turbidite system (western North America) and distal basin-floor deposits of the Ordovician Cloridorme Formation (eastern North America). Stacking patterns derived from field-based observations are then evaluated using a first-order Markov chain (e.g., Krumbein, 1968; Carr, 1982; Davis, 1986) to quantitatively determine whether there is any systematic stacking arrangement, and if so, identify possible physical drivers responsible for the pattern. This statistical technique has previously been applied to describe and interpret deep-marine matrix-rich strata at either the individual bed scale (e.g., Lowe and Guy, 2000; Haughton et al., 2003; Davis et al., 2009) or the stratal element scale (e.g., Kane and Pontén, 2012; Terlaky and Arnott, 2016); however, in this study, Markov chain analysis is carried out to assess the stacking patterns of beds and stratal units of matrix-rich and associated matrix-poor lithologies, and this study therein attempts to address the gap in understanding that remains between the deposition of an individual bed and that of a stratal element.

The Neoproterozoic (740–570 Ma) Windermere Supergroup (WSG) is an unconformity-bounded succession of mostly metasedimentary rocks deposited on the passive continental margin of Laurentia (ancestral North America) during the breakup of the Rodinia supercontinent (Stewart, 1972; Ross and Arnott, 2007; McMechan, 2015). The outcrop belt of the WSG stretches from northwestern Mexico through the western United States and Canada and northward to the Yukon-Alaska border (Ross and Arnott, 2007) (Fig. 1A). In the southern Canadian Cordillera, the basal rift sequence includes a few-kilometers-thick assemblage of intercalated glacial diamictites and mafic volcanic rocks of the Toby and Irene Formations, respectively (Aalto, 1971; Ross et al., 1995; Warren, 1997). This is then overlain by a 5–7-km-thick post-rift, upward-shoaling succession composed of deep-marine basin floor (Kaza Group) followed by continental slope (Isaac Formation) to upper continental slope and shelf (Cunningham and Yankee Belle Formations, respectively) deposits (Campbell et al., 1973).

Deep-marine strata of the WSG basin (Fig. 1B), termed the Windermere turbidite system, are well exposed throughout the southern Canadian Cordillera (Fig. 1A). At the Castle Creek study area (Fig. 1D), vertically dipping, recently deglaciated, vegetation-free strata of the Upper Kaza Group and Isaac Formation form an outcrop belt that is ~2.6 km thick and 8 km wide (Terlaky et al., 2016; Navarro and Arnott, 2020). Here, field observations of beds and bedsets can be carried out over distances of tens of meters vertically, but more significantly, hundreds of meters along the strike of the outcrop. However, due to the two-dimensionality of the outcrop, paleoflow data are limited to a small number of well-exposed three-dimensional current ripples and dune cross-stratified sets and a small number of flute marks, which collectively indicate a general transport direction ranging from west-northwest to northeast (see Ross and Arnott, 2007; Schwarz and Arnott, 2007; Khan and Arnott, 2011; Navarro and Arnott, 2020) and therefore subparallel to the regional paleoflow direction (e.g., Mountjoy and Aitken, 1963).

Furthermore, rocks in the study area have undergone lowgrade (greenschist facies) metamorphism, resulting in the alteration of primary clay matrix to muscovite and chlorite and strain-induced structures like bulging recrystallization and sub-grain rotation along some grain boundaries. Despite these changes, primary sedimentary structures and most sedimentary textures, namely grain size, grain shape, and sorting, are well preserved, enabling these rocks to be described and classified as sedimentary rocks (Ross and Arnott, 2007; Ningthoujam et al., 2022).

The Middle Ordovician (460–445 Ma) Cloridorme Formation (CF) comprises the Quebec segment of a succession of turbidite systems that were deposited along an elongate foreland basin located between cratonic North America to the northwest and the developing Taconic orogen and its associated volcanic arcs to the southeast (St. Julien and Hubert, 1975; Hiscott et al., 1986). The basin extended from Newfoundland (eastern Canada) to Tennessee (eastern United States) and was as much as 150 km wide (Hiscott, 1984). The CF is exposed semi-continuously for ~150 km along the north shore of the Gaspé Peninsula, Quebec, Canada (Fig. 1A). The CF is at least 4 km thick and, based on lithofacies and extensive marker beds (megaturbidites), is divided, from oldest to youngest, into five members: Saint Hélier member, composed of distal basin-floor mudstones; Pointe-à-la-Frégate and Petite-Vallée members, consisting of sandstone-rich lobe and mudstone-rich inter-lobe deposits; Mont-Saint-Pierre member, made up of distal basin-floor mudstones; and Marsoui member, made up of sand lobes and interlobe deposits (Hiscott et al., 1986) (Fig. 1C). This study focuses on strata of the Petite-Vallée member located near the village of Petite Vallée, Quebec (Fig. 1E). Here, strata are nearly vertically dipping (average 78°) and east-to-west (271°) striking and exposed along an intertidal wave-cut platform that allowed for field observations of beds and bedsets to be carried out over tens of meters vertically but also tens to hundreds of meters along the strike of the outcrop. In the study area, average paleocurrent based on base-of-bed tool and flute marks and troughs of three-dimensional current ripples is 272° (n = 102) and therefore sub-parallel to the strike of the outcrop, which is consistent with previous studies in other parts of the CF (Parkash and Middleton, 1970; Pickering and Hiscott, 1985; Awadallah, 2002). Rocks in the study area have experienced very lowgrade (zeolite to prehnite-pumpellyite facies) metamorphism associated with the Taconic orogeny and resulted in the alteration of primary clay minerals to illite and chlorite (Jiang and Peacor, 1994). Despite this change, primary sedimentary structures and textures are well preserved, allowing these strata to be described and classified using conventional sedimentary terminology (see Ningthoujam et al., 2022).

In slope and proximal basin-floor deposits of the WSG and distal basin-floor deposits of the CF, four main lithofacies linked with matrix-rich and associated matrix-poor strata were recognized (Terlaky and Arnott, 2014; Angus et al., 2019; Ningthoujam et al., 2022). Based on matrix content, the four lithofacies are: matrix-poor sandstone (MPS) with 0%–20% matrix content; muddy sandstone (MS) with 20%–50% matrix content; bipartite facies (BF) with a basal sandy (30%–60% matrix) part overlain sharply by a planar- to irregular-based muddier upper part (40%–80% matrix); and sandy mudstone (SM) with 50%–90% matrix content (Fig. 2). All four facies are structureless (i.e., lack tractional sedimentary structures), massive or coarsetail graded, and composed of rounded to angular, equant to elongate framework grains that range from very fine to very coarse sand. Mud and lesser very fine- to medium-grained sandstone clasts are also common. Clasts are typically elongate, sub-angular to sub-rounded, and oriented with their maximum projection plane sub-parallel to bedding. Where present, these strata are then capped by lamina to thin beds of traction-structured sandstone overlain by silty mudstone. Tracing of individual beds along the strike of the outcrop in the WSG (Terlaky and Arnott, 2014; Angus et al., 2019; Wearmouth, 2018; Ningthoujam et al., 2022) and CF (Ningthoujam et al., 2022) reveals a systematic arrangement of lithofacies, namely, MPS to MS to BF and then to SM, each overlain, except where eroded, by a laminated to thinbedded, traction-structured sandstone and silty mudstone (Fig. 2). More specifically, thick-bedded MPS gradually transforms to medium- to thick-bedded MS, which then develops into medium-bedded two-part BF consisting of a basal muddy sandstone overlain sharply by an upper sandy mudstone commonly containing dispersed mud clasts. Further along the transect, the basal part of BF progressively thins and then pinches out, whereas the upper part thickens but at a diminishing rate. Further laterally, following the termination of the basal sandy part of BF, the now exclusively SM thins to a pinch-out (Fig. 2). Importantly, these changes occurred over a horizontal distance of a few hundreds of meters and were interpreted to reflect a downflow (proximal to distal) depositional continuum, which in the WSG was interpreted to be oriented at a high angle to the main flow (Terlaky and Arnott, 2014; Wearmouth, 2018; Angus et al., 2019; Ningthoujam et al., 2022) but more or less parallel to the main flow in the CF (Ningthoujam et al., 2022). Significantly, regardless of the orientation of the depositional transect (parallel or at high angle) to the main flow or of depositional setting (continental slope versus proximal or distal basin floor), the transitions from MS to BF, BF to SM, and SM to its pinch-out were reported to be of proportionate length (Ningthoujam et al., 2022) and, accordingly, suggest similarity in depositional process.

Based on field observations, matrix-rich and associated matrix-poor strata stack exclusively to form stratal elements several meters to tens of meters thick that typically underlie sand-rich stratal elements such as slope channels, distributary channels, and terminal splays in the continental slope to proximal basin floor settings of the WSG (Terlaky and Arnott, 2014; Terlaky et al., 2016; Angus et al., 2019; Ningthoujam et al., 2022) (e.g., Fig. 3) and sandy terminal splays in the distal basin floor setting of the CF (Ningthoujam et al., 2022) (e.g., Fig. 4). This stratigraphic relationship was attributed to avulsion, and the stratal elements were identified as “avulsion splays” (sensuTerlaky and Arnott, 2014). Avulsion formed high-energy, momentum-driven turbulent wall jets that eroded the mud-rich inter-channel area and charged the flow with fine-grained sediment, namely very fine sand to clay, in addition to low-density mud clasts (Arnott, 2007a, 2007b; Terlaky and Arnott, 2014; Angus et al., 2019; Navarro and Arnott, 2020; Ningthoujam et al., 2022). The local incorporation of easily suspended fine-grained sediment resulted in an abrupt increase in local effective viscosity and a dramatic change in the characteristics and reduction in the diffusion of turbulence required to maintain the particle suspension. This caused the sediment suspension to rapidly exceed transport capacity and collapse, forming a negligibly sheared suspension that deposited a systematic along-flow depositional continuum consisting of MPS to MS to BF and then SM over distances of hundreds of meters. Following the emplacement of the depositional continuum from the collapsed suspension, a remnant low-density turbidity current reworked the top of the deposit and formed a thin-bedded, very fine-grained, upper-division turbidite capped by a silty mudstone deposited by a gelled near-bed fluid-mud layer. In continental slope and proximal basin-floor deposits of the WSG, the depositional transect developed along the margins of the high-energy wall jet, whereas in the distal basin-floor setting of the CF, deposition took place in the same direction as the main flow and formed a mouth bar that marks the downflow terminus of the flow (see discussion and figure 17 in Ningthoujam et al., 2022).

To better elucidate the stacking pattern of matrix-rich and associated matrix-poor strata, a combination of field-based observations and a statistical analysis was used. Field observations include stratal thickness, bedding contacts, grain size, and sedimentary structures in each bed. Hand samples of matrix-rich and matrix-poor lithologies were collected in the field and then analyzed using an optical petrographic microscope. Individual thin sections were point counted (≥300 points per sample) to accurately describe textural characteristics, and the longest dimension of each framework grain (grain size between 0.0625 mm and 6 mm) was measured to determine grain-size distribution in each sample. Notably, all the petrographic work was conducted by a single operator to minimize operator bias and standardize potential error.

Based on earlier stratigraphic and depositional-framework studies at Castle Creek (WSG), select intervals containing matrix-rich and associated matrix-poor lithologies were chosen for this study. In proximal basin-floor deposits of the Upper Kaza Group, a 100-m-thick and 1000-m-wide stratigraphic interval made up of a succession of terminal splays, distributary channel fill, debrites, thin-bedded turbidites, and matrix-rich and associated matrix-poor strata mapped previously by Terlaky and Arnott (2014, 2016), Terlaky et al. (2016), Angus et al. (2019), and Ningthoujam et al. (2022) was selected. In this section, 61 stratigraphic logs (5–50 m long) were measured over a lateral distance of ~900 m targeting matrix-rich and associated matrix-poor strata (Fig. 3A). Similarly, in slope deposits of the stratigraphically higher Isaac Formation, two stratigraphic intervals (~80 m thick and 400 m wide, and ~120 m thick and 400 m wide) comprising sandy channel fills, thin-bedded turbidites, and matrix-rich and associated matrix-poor strata mapped previously by Arnott (2007a, 2007b), O’Byrne et al. (2007), Angus et al. (2019), Arnott et al. (2021), and Ningthoujam et al. (2022) were selected. Here, a total of 32 logs (3.5–30 m long) were measured in two transects targeting matrix-rich and associated matrix-poor strata (Figs. 3B and 3C). In distal basin-floor deposits of the Petite Vallée member (CF), an ~300-m-thick and 1700-m-wide stratigraphic interval comprising terminal splays, thin-bedded turbidites, and matrix-rich and associated matrix-poor strata mapped previously by Enos (1969), Beeden (1983), Hiscott et al. (1986), and Ningthoujam et al. (2022) was selected. Here, 15 bed-by-bed stratigraphic logs (18–147 m long) were measured (Fig. 4). Also, in both the WSG and the CF, a few logs were extended upward to capture details of the overlying stratal element (Figs. 3 and 4). Additionally, hand samples of matrix-rich and associated matrix-poor lithologies were collected at both study sites, and a total of 253 standard (30 µm thick) thin sections were prepared for conventional petrographic analysis and point counting.

A first-order Markov chain analysis (e.g., Terlaky and Arnott, 2016) was conducted at bed and stratal-unit scales to quantitatively assess how matrix-rich and associated matrix-poor strata stack to build up a stratal element. This statistical technique describes sequential changes in a system composed of discrete, finite states (Krumbein, 1968; Schwarzacher, 1975). The change from one state to the next is considered a transition and is assumed to depend only on the state that precedes it, an assumption termed a first-order Markov property. Also, the system is considered to change randomly from one state to another based on the transition probabilities between states. Furthermore, the stratigraphic sections used in this statistical analysis are assumed to be deposited under static paleoenvironmental conditions with no changes due to external forcing (e.g., climate, tectonism, sediment supply) (e.g., Lehrmann and Goldhammer, 1999; Terlaky and Arnott, 2016). Lastly, the vertical transition of states (i.e., from one bed to the next bed, or one stratal unit to the next stratal unit) observed in the stratigraphic sections assumes no erosion and that all states are fully preserved. In the WSG, the input data were obtained from a composite stratigraphic log consisting of all measured matrix-rich and associated matrix-poor strata in the Castle Creek study area. Here, individual logs with the most well-preserved stratigraphy were chosen from each section to build a single continuous 137-m-long stratigraphic section (Fig. 5). In the CF, on the other hand, data were derived from a single, vertically continuous stratigraphic log (286 m long) (Fig. 6). The observed vertical stacking of these strata was then used to construct the transition matrices, which then were tested against the expected transition matrices (zero-order Markov chain, or the null hypothesis) to evaluate statistical significance.

In this study, a stratigraphic hierarchy is recognized. At the most basic level, a “bed” consists of one of the four lithofacies recognized by Ningthoujam et al. (2022), namely, matrix-poor sandstone (MPS), muddy sandstone (MS), bipartite facies (BF), and sandy mudstone (SM), which, except where eroded, is then capped by a <1-cm-to 10-cm-thick layer of traction-structured sandstone overlain by a silty mudstone. Each bed is then named by the principal lithofacies type (e.g., a bed with BF is called bipartite facies bed) (Fig. 7). Beds then stack to form “units”, each of which consists of two or more beds of similar lithofacies type or, less commonly, a single bed that is overlain and underlain by lithologically dissimilar strata (Fig. 7). Additionally, single-bed stratal units occur as single layers intercalated with thin-bedded upper-division turbidites (Fig. 6). Like individual beds, stratal units are named according to the principal lithofacies type (e.g., a stratal unit comprising BF is called a bipartite facies stratal unit). Stratal units then stack to form a stratal “element”, which, being composed of either exclusively matrix-rich and associated matrix-poor lithologies or matrix-rich strata intercalated with subordinate thin-bedded turbidite units, is interpreted as an “avulsion splay” (sensuTerlaky and Arnott, 2014) (Figs. 5–7).

Based on field observations in both study areas, beds of similar litho-facies, but of variable thickness, typically stack (70% of all measured beds in WSG and 65% in CF) to form multi-bed stratal units consisting of two to nine beds (e.g., Figs. 5 and 6). Less commonly (30% of all measured beds in WSG and 35% in CF), the succession consists of single-bed stratal units, which in some places stack but are differentiated from turbidites or other single beds based on unique textural and dimensional attributes (e.g., Figs. 5 and 6). Stratal-unit boundaries are abrupt and typically marked by a change in lithofacies (Figs. 5–7). Beds within a multi-bed stratal unit undergo a transformation from one facies to the next in the lateral depositional continuum, for example, MS to BF, or BF to SM, at about the same lateral position (Fig. 8). Additionally, in the WSG, at any position along a multi-bed stratal unit, the D90 (coarse fraction, 90th percentile of particle size) and D50 (median particle size) of the grain-size distribution subtly and progressively decrease upward (Fig. 9). In the CF, on the other hand, grain-size distribution shows a slight but nonetheless abrupt decrease in the topmost bed in each multi-bed stratal unit (Fig. 10). In both study areas, matrix content changes negligibly upward within a multi-bed stratal unit (e.g., Figs. 9 and 10).

A total of 333 and 304 bed transitions were observed in CF and WSG, respectively. Results of the observed and expected transitions of beds in both study areas are summarized in the transition matrix and state diagrams (Fig. 11). In WSG, MPS beds comprise 36% of all beds (i.e., states) observed and are most frequently overlain by similar MPS beds (48% of observed transitions). Less commonly, MPS beds are overlain by beds of SM (27%), MS (20%), or BF (5%). MS beds make up 15% of all the observed beds and most commonly are overlain by other MS beds (34%), and less commonly by beds of SM (28%), MPS (25%), or BF (13%). BF beds constitute 17% of all the beds observed and are most frequently overlain by other BF beds (48%) followed by beds of MPS (25%), SM (21%), or MS (6%). Lastly, SM beds comprise 32% of all the observed beds and are usually overlain by other SM beds (44%), and less commonly by beds of MPS (34%), BF (16%), or MS (6%).

In CF, MPS beds make up 12% of all observed beds and most commonly are overlain by beds of SM (32%) and BF (32%) and less commonly by beds of MPS (29%) or MS (7%). MS beds comprise 16% of all the observed beds and are most frequently overlain by BF beds (48%). Less frequently, MS beds are overlain by other MS beds (25%) and beds of SM (17%) or MPS (10%). BF beds constitute 41% of all observed beds and most commonly are overlain by other BF beds (50%) followed by beds of SM (26%), MS (15%), or MPS (9%). Lastly, SM beds make up 31% of all the observed beds and most frequently are overlain by other SM beds (44%), and less frequently by beds of BF (29%), MS (15%), or MPS (12%).

In the statistical analysis, the degrees of freedom are (n – 1)2 = (4 – 1)2 = 9, where n is the number of discrete states (i.e., MPS, MS, BF, and SM beds), and the minimum chi-squared (χ2) value to reject the null hypothesis (at 5% significance level) is 16.92. The values in the χ2 matrix are calculated using χ2 = (observed – expected)2 / expected. The respective sums of the χ2 matrices are 76 for WSG and 33 for CF (Figs. 11C and 11G). Accordingly, the null hypothesis can be rejected for both study areas and the observed vertical stacking of beds is non-random.

More specifically, comparing the expected and observed transition matrices in both study areas shows that beds of similar lithofacies type typically succeed one another, whereas all other transitions are either less than expected, random, or very close to random (Figs. 11C, 11D, 11G, and 11H). Importantly, in CF, for example, a statistically significant transition does not necessarily correlate with a high or low transition percentage. For example, the transitions from MPS bed to BF bed and MPS bed to SM bed are both 32%, whereas MPS bed to another MPS bed is 29%. Despite the first two transitions being more frequent than the third, the first two are statistically very close to random or random, whereas the third is statistically more likely than random (Fig. 11G). This shows that the most or the least likely transitions in terms of absolute number of occurrences are not necessarily statistically significant and in fact may be only the result of random chance.

In both study areas, stratal units of matrix-rich and associated matrix-poor lithologies stack exclusively to form 0.4–33.2-m-thick (average 2.7 m) avulsion splays (Figs. 5 and 6). These avulsion splays are equivalent in hierarchy to “lobe elements” of Prélat et al. (2009) and in the WSG abruptly underlie sandrich stratal elements such as slope channels (e.g., Arnott, 2007a, 2007b; Schwarz and Arnott, 2007; Angus et al., 2019; Ningthoujam et al., 2022) (e.g., Fig. 3), distributary channels and terminal splays in the channel-lobe transition zone (Navarro and Arnott, 2020), and proximal basin floor (Terlaky and Arnott, 2014; Terlaky et al., 2016; Angus et al., 2019; Ningthoujam et al., 2022) (e.g., Fig. 3). In distal basin-floor strata of the Petite Vallée study area, the avulsion splays typically underlie terminal splays (see Ningthoujam et al., 2022) (e.g., Fig. 4).

A total of 203 and 184 stratal unit transitions were observed in the CF and WSG, respectively. Results of the observed and expected transitions of stratal units in both study areas are summarized in the transition matrix and state diagrams (Fig. 12). In the WSG, MPS stratal units comprise 34% of all stratal units observed and are most frequently overlain by stratal units of SM (52% of observed transitions). Less commonly, MPS stratal units are overlain by stratal units of MS (30%), BF (14%), or MPS (4%). MS stratal units make up 17% of all observed stratal units and most commonly are overlain by stratal units of SM (50%) and less commonly by MPS (38%) or BF (12%). Furthermore, MS stratal units were not observed to overlie one another. Stratal units of BF constitute 14% of all stratal units observed and are most frequently overlain by stratal units of SM (46%) or MPS (46%), and less commonly by MS (8%). Additionally, BF stratal units were not observed to overlie one another. Lastly, SM stratal units comprise 35% of all observed stratal units and most commonly are overlain by stratal units of MPS (59%) and less commonly by BF (21%), MS (17%), or SM (3%).

In the CF, MPS stratal units make up 12% of all observed stratal units and most commonly are overlain by stratal units of SM (54%) and less commonly by BF (38%) or MS (8%). Additionally, MPS stratal units were not observed to overlie one another. Stratal units of MS comprise 21% of all observed stratal units and are most frequently overlain by BF stratal units (62%). Less frequently, MS stratal units are overlain by stratal units of SM (24%), MPS (7%), or MS (7%). BF stratal units constitute 36% of all observed stratal units and most commonly are overlain by stratal units of SM (47%) followed by MS (26%), MPS (15%), or BF (12%). Lastly, SM stratal units make up 31% of all the observed stratal units and are most frequently overlain by stratal units of BF (45%) and less frequently by MS (28%), MPS (16%), or SM (11%).

In the statistical analysis, the degrees of freedom are (n – 1)2 = (4–1)2 = 9, where n is the number of discrete states (i.e., MPS, MS, BF, and SM stratal units) and the minimum chi-squared (χ2) value to reject the null hypothesis (at 5% significance level) is 16.92. The respective sums of the χ2 matrices are 76 for WSG and 53 for CF (Figs. 12C and 12G). Therefore, the null hypothesis can be rejected for both study areas and the observed vertical stacking of stratal units is non-random.

Notably, comparison of the expected and observed transition matrices in both study areas shows that similar stratal units are less likely to overlie themselves than expected (Figs. 12C, 12D, 12G, and 12H). With the exception of MPS to MS, MPS to SM, and SM to MPS in the WSG, and MPS to SM, MS to BF, and BF to SM in the CF, all other stratal unit transitions are either random or very close to random (Figs. 12C, 12D, 12G, and 12H).

Based on field observations and verified statistically, beds of similar litho-facies, but commonly of different thickness, typically succeed one another vertically and form multi-bed (two to nine beds) stratal units; all other transitions are either less than expected, random, or very close to random (Fig. 11). Within a discrete multi-bed stratal unit, successive beds transition from one facies to the next at about the same lateral position (Fig. 8). Also, grain size shows a subtle but nevertheless consistent upward decrease; matrix content changes negligibly upward (Figs. 9 and 10). Collectively, these observations suggest that beds within a multi-bed stratal unit are more or less texturally self-similar and therefore deposited from flows with similar hydraulic and compositional character, albeit exhibiting a minor but nevertheless progressive decrease in flow energy. Superimposed on this temporal regularity is the deposition of the traction-structured turbidite overlain by mudstone that caps each bed. Each cap signals a dramatic change in depositional conditions, which then revert to conditions that resembled those that deposited the underlying bed, albeit commonly with a slightly finer-grained sediment supply. This, then, raises the obvious question as to the origin of such systematic sedimentation and whether deposition of a multi-bed stratal unit is a consequence of multiple flow events or one flow event with multiple pulses.

With respect to multiple flow events, it has been shown experimentally that oversized turbidity currents traversing a partially confined sinuous channel equilibrate to the channel’s cross-sectional area by flow stripping and overspill in channel bends, which ultimately reduces size variability in the flows (Straub et al., 2008). More recently, Kelly et al. (2019) described a process termed flow tuning where saline currents passing through a partially confined straight channel adjusted their flow size either by overspilling (oversized flows) or by inflating through ambient fluid entrainment (undersized flows) to achieve an equilibrium flow condition in which overbank losses are balanced with ambient fluid entrainment. Intriguingly, the self-similarity of beds within a multi-bed stratal unit could be consistent with successive flows that had undergone flow tuning; however, variability in the thickness of individual beds within a multi-bed stratal unit (e.g., Figs. 5 and 6) and the subtle but nonetheless consistent upward decrease of grain size observed in multi-bed stratal units makes this problematic (Figs. 9 and 10).

A second possibility is that self-similar beds were deposited from multiple pulses of a single turbidity current. Pulsing in turbidity currents has been recognized in direct monitoring of natural flows (e.g., Shepard et al., 1975; Best et al., 2005; Menczel and Kostaschuk, 2013) and studied experimentally and simulated numerically (e.g., Dai, 2008; Ho et al., 2018a, 2018b; Kostaschuk et al., 2018); however, relatively few studies have described evidence of this in the ancient sedimentary record. In Lowe (1982) and Larue and Provine (1988), stacks of self-similar turbidites lacking evidence of intervening erosion were suggested to have been formed by pulsations or surges in velocity during a single flow event. Additionally, alternating arenite and pelite turbidite facies, equivalent to stacked Bouma Tcd/e (upper division turbidite beds) divisions, were reported from Oligocene and Carboniferous rocks in eastern France by Duringer et al. (1991) and interpreted to have been deposited by a single, continuous muddy turbidity current with episodic higher-energy pulses forming the Tc divisions. Also, in the Precambrian Kongsfjord Formation in Norway, Pickering (1979) interpreted ~2–5-m-thick bedsets, each comprising at least five beds and composed mostly of upper-division turbidites and exhibiting a progressive upward fining and thinning, to record pulses of traction-transport deposition separated by lower-energy suspension deposition during a single continuously waning flow event. Most recently, Cunningham and Arnott (2021) described thin-bedded, upper-division turbidites in continental-slope levee deposits in the WSG that form bedsets of about two to ten self-similar beds, which they attributed to deposition from multiple episodic pulses of a single channelized flow event.

Although it seems clear from these outcrop examples that stacked bed self-similarity is suggestive of turbidity current pulsing, the cause for the pulsing still remains poorly understood. To date, a variety of mechanisms have been proposed to explain pulsing, including: high swell waves and high tides generating pulsating turbidity currents (Shepard et al., 1975); retrogressive slumping that evolved downslope into multi-pulsed turbidity currents (Piper et al., 1999; Canals et al., 2004; Ho et al., 2018a); multi-pulsed flows formed downflow of channel confluences by discrete flows triggered simultaneously (e.g., by an earthquake) in different upflow channels such as tributaries in a canyon system (Nakajima and Kanai, 2000; Ho et al., 2018a); velocity pulsing caused by hydrodynamic instabilities in plunging hyperpycnal flows (Kostaschuk et al., 2018); hypo- and mesopycnal flows generating multiple secondary turbidity currents through settling-driven Rayleigh-Taylor convection (Hoyal et al., 1999; Parsons et al., 2001; Davarpanah Jazi and Wells, 2016, 2020). Additionally, in some confined settings, like in the CF, flow reflection has been suggested to cause pulsing in turbidity currents (e.g., Haughton, 1994; Ho et al., 2018a). Most of these mechanisms, except for retrogressive slope failures, are expected to have short delay times (few minutes to hours) between pulses, which would preclude the emplacement of the traction-structured sandstone and mud cap over individual beds in a bedset. In a multi-pulsed flow initiated by retrogressive slumping, the time between successive slope failures can be sufficient to allow for the deposition and preservation of traction-structured sandstone and mud cap, but it is highly unlikely to result in consistent self-similarity of successive pulses of sediment transport and deposition. Additionally, the subtle but nevertheless consistent upward-fining trend observed in all multi-bed stratal units is difficult to reconcile with any of these mechanisms (Figs. 9 and 10).

An alternative mechanism builds on recent direct monitoring of turbidity currents in the Congo Canyon by Azpiroz-Zabala et al. (2017). Unlike typical laboratory-generated surge-type turbidity currents, where the body is faster than the head and continuously feeds the head with sediment (e.g., Middleton, 1966), Azpiroz-Zabala et al. reported that Congo Canyon turbidity currents consist of a dense, fast-moving and erosive head, termed a “frontal cell”, and a trailing body and tail (Fig. 13A). Here, the flow is sustained by the frontal cell through continuous erosion and incorporation of new sediment from the sea-floor while simultaneously shedding sediment-laden fluid into the trailing part of the flow. The difference in density and hence velocity between the frontal cell and the rest of the flow (i.e., the body and tail) allows the frontal cell to outrun the trailing flow and stretch the flow downstream. This process also allows turbidity currents to last as long as 10 days compared to a few minutes to hours in previous flow measurements in other modern systems (e.g., Xu et al., 2004; Khripounoff et al., 2012; Liu et al., 2012; Hughes Clarke, 2016). Considering that the flow monitoring by Azpiroz-Zabala et al. was limited to the submarine canyon and that at least some turbidity currents are capable of reaching further downslope and even onto the basin floor, it is possible that as the turbidity current continues to be stretched, the dense frontal cell eventually becomes detached from the rest of the flow and forms a separate turbidity current (Fig. 13B). With the detachment of the head, the remaining (main) turbidity current is cut off from its major source of sediment, which then causes it to be diluted and slow. Additionally, elevated stagnation pressure at the new leading edge of the main flow causes it to move more slowly than the trailing part, which then feeds sediment forward and into a newly developing head (Middleton, 1966; Kneller and Buckee, 2000; Cartigny et al., 2013). With time, increased density in the head causes it to accelerate and erode the seafloor, which then feeds new sediment into the head and also sheds into the trailing body. Farther downflow, the head evolves into a new frontal cell (Fig. 13C) that eventually detaches from the main flow and becomes a second discrete turbidity current (Fig. 13D). With time and distance downslope, this process is possibly repeated several times, creating multiple surges from a single flow event. Accordingly, each bed in a multi-bed stratal unit would correlate with a single surge (Fig. 13E), and the progressive decrease in grain size reflects the gradual longterm waning of the parent turbidity current. Additionally, single-bed stratal units that separate multi-bed stratal units or occur as a single layer in a succession of thin-bedded turbidite units (Figs. 5 and 6) may indicate flows where only the initial surge (i.e., first detached head) reached that downslope position while the rest of flow terminated some distance further upslope.

Importantly, the delay time between successive surges must have been sufficient to allow for the emplacement of traction-structured, very fine-grained sandstone (0.2–3 cm thick) and overlying mud cap. Similar low-frequency pulses have also been suggested by Cunningham and Arnott (2021) to allow for the deposition and semi-consolidation of silty and muddy Td/e divisions overlying self-similar sand-rich Tc divisions in continental-slope levee deposits of the Isaac Formation (WSG). Potentially more significant is the mud cap, which typically is massive and 3–5 cm thick. Here the rapid accumulation of fine-grained sediment is attributed to very high near-bed sediment concentration and extensive particle-particle interaction that promoted mud flocculation and aggregation (McAnally et al., 2007; see also discussion in Cunningham and Arnott, 2021). Particle aggregation may also have been enhanced by mud particle adhesion to suspended microbial extra-cellular polymeric substances (EPS) (Bennett et al., 1992). Given the absence of bioturbation in WSG and CF strata, it is possible that the mud-rich seabed was anomalously enriched in microbes and associated EPS, which, when eroded during avulsion, resulted in rapid and extensive particle aggregation in the near-bed region. Eventually, the concentration of flocculated and aggregated particles led to gelling and en masse deposition of the mud layer. Notably, the contact between a mud cap and the overlying bed in a multi-bed stratal unit is typically sharp and planar. This indicates that the mud caps had sufficient strength and load-bearing capacity to resist shearing and loading during deposition of the overlying bed. If strength were to have developed through self-weight consolidation (i.e., normal consolidation; Skempton, 1969), the time frame would have greatly exceeded the predicted delay time between surges associated with a single flow event. However, based on observations of modern deep-marine sediments in offshore Angola and Nigeria, the Nova Scotian Slope, and in the Pacific Ocean (Baltzer et al., 1994; Meadows and Meadows, 1994; Ehlers et al., 2005; Kuo and Bolton, 2013), a zone of exceptionally high strength, albeit with high water content, and termed the “crust” (sensuKuo and Bolton, 2008) occurs a few tens of centimeters to about a meter below the seafloor with strength that is an order of magnitude higher than from normal consolidation. Although development of the crust is still poorly understood, microbially mediated strength enhancement of marine sediments through production of EPS has been proposed in several studies (e.g., Dade et al., 1990; Tolhurst et al., 2002; Bhaskar and Bhosle, 2005). Additionally, Parkes et al. (2000) showed that the concentration of bacteria in the uppermost meter of the marine sediment pile can exceed 109 cells per cubic centimeter of sediment and rival that in the terrestrial realm. Using similar concentration of microbes (i.e., 109 cells per cubic centimeter of sediment) as measured in modern subsurface marine sediments, Dade et al. (1996) experimentally showed that a glue-like exopolymer produced by marine bacteria (Alteromonas atlantica) enhanced the yield stress of a kaolinite clay and seawater mixture by 60% in ~0.5–9 days. Collectively, it is postulated that a lag of possibly several days between successive surges, combined with biophysically enhanced mud particle aggregation in the near-bed region, led to gelling, and en masse deposition of the mud layer, as well as microbe-induced post-depositional strength enhancement of the mud layer, allowed for the deposition and preservation of mud caps associated with successive surges during the same flow event.

Results of the Markov chain analysis of matrix-rich and associated matrix-poor lithologies at the stratal unit scale show that similar lithologies that make up part of the depositional continuum in each of their respective stratal units stack less commonly than expected in a random distribution in both study areas, whereas stratal unit transitions of MPS to SM, SM to MPS, and MPS to MS in the WSG and MS to BF, BF to SM, and MPS to SM in the CF are observed significantly more frequently than expected in a random distribution (Figs. 12C, 12D, 12G, and 12H). This indicates that the observed vertical stacking in both study areas is more heterogeneous than predicted by a random distribution of stratal units, and the fact that this statistical analysis (see Methods and Terminology section above) assumes no external forcing (e.g., variations in sediment supply, changes in relative sea level, tectonic variability) but also the commonality of these distinctive strata in the deep-water stratigraphy of both the WSG and CF imply an internal forcing.

In WSG strata, the stack of the thickest and most proximal lithofacies in a single stratal unit (i.e., MPS stratal unit) is preferentially overlain by a stack of the thinnest and most distal lithofacies in the overlying stratal unit (i.e., SM stratal unit), or vice versa (Figs. 12C and 12D). Considering that these strata were deposited on opposite ends of the depositional continuum, the stacking pattern indicates an abrupt and long-distance lateral displacement of facies. This is the stratigraphic manifestation of compensational stacking, which here resulted from the flows being diverted around topographic highs formed by deposition of the previous stratal unit and directed into adjacent topographic lows (Terlaky and Arnott, 2016) (Fig. 14A). Flow diversion was followed by a temporary stabilization of flow and sedimentation conditions and deposition of the next stratal unit. MPS stratal units are also preferentially overlain by slightly more distal MS stratal units, however the opposite transition (i.e., MS to MPS) is shown not to be statistically significant (χ2 value of 0.1) (Figs. 12C and 12D). The MPS-overlain-by-MS transition indicates a minor backstep of successive stratal units that might reflect flow deceleration against the positive relief created by deposition of the underlying stratal unit (Pratson et al., 2000; Ferguson et al., 2020) (Fig. 14A). In the CF strata, on the other hand, preferential stacking of stratal units including MS to BF, BF to SM, and MPS to SM are observed (Figs. 12G and 12H). Notably, however, the inverse of these transitions shows no statistical significance, suggesting a strong preference for systematic backstepping at the stratal unit scale. As noted earlier, matrix-rich and associated matrix-poor strata in the CF form mouth-bar deposits at the downflow terminus of the flow (see Ningthoujam et al., 2022); therefore, this systematic back-stepping of stratal units most likely represents accretion on the upflow side of the mouth bar (Fig. 14B).

Overall, results of the statistical analysis suggest a preferred growth pattern of stratal units consisting of initiation followed by gradual back-stepping and then side-stepping into an adjacent topographic low (WSG) or consistent back-stepping (CF). The stacking pattern can either occur once or repeat multiple times building up several-meters-thick stratal elements termed avulsion splays (Terlaky and Arnott, 2014). A similar growth pattern is commonly reported from submarine fans and at a variety of scales, including beds and/or bedsets (e.g., Deptuck et al., 2008; Jobe et al., 2017), lobe elements (e.g., Deptuck et al., 2008; Jegou et al., 2008), and lobes (e.g., Twichell et al., 1991; Vittori et al., 2000; Fildani et al., 2018). The cause for this systematic growth pattern is still a topic of debate, and studies have associated it with either allogenic controls, like sea level and climate (Hodgson et al., 2006; Marsset et al., 2009; Romans et al., 2009; Jobe et al., 2017), or internal autogenic controls, like the interaction between flows and the evolving topography (Gervais et al., 2006; Jegou et al., 2008; Prélat et al., 2009). Although allogenic forcing may exert some influence on the depositional pattern of the packages described here, Cantelli et al. (2011) and Fernandez et al. (2014) showed experimentally that similar depositional patterns (i.e., forward-stepping, back-stepping, and then side-stepping of the depocenter) can be generated with no change in allogenic conditions. Additionally, more recent numerical modeling (Burgess et al., 2019) and physical experiments (Ferguson et al., 2020) showed that in unconfined settings like the basin floor, even when external factors change, being sediment supply in these studies, the effect of depositional topography (i.e., autogenic forcing) on the evolution and orientation of the flows dominated and generated depositional patterns that were largely similar to those in earlier studies with constant sediment supply.

Even though results of the statistical analysis show that the stacking pattern of lithofacies stratal units making up an avulsion splay are similar to many of the bed- and/or bedset- to lobe-scale stacking patterns recognized in both modern and ancient deep-marine systems, the composite strata are unique in their high mud content, which then requires an abrupt and temporarily sustained supply of mud. As interpreted above, and previously by Terlaky and Arnott (2014), Angus et al. (2019), and Ningthoujam et al. (2022), the abrupt development of matrix-rich strata is caused by an upflow avulsion that activates the local transport system. During avulsion, the flow loses confinement and forms a wall jet that locally scours the muddy inter-channel area, thereby charging the flow with mud (see Ningthoujam et al., 2022). The scour locally steepens the seafloor gradient, forming an upflow-migrating knickpoint that locally increases the erosive potential of successive flows (e.g., Fildani et al., 2013; Carvajal et al., 2017). Upflow-migrating knickpoints have become increasingly recognized in studies of the modern seafloor and shallow subsurface, especially those using high-resolution seafloor imaging techniques (e.g., Fildani et al., 2006, 2013; Armitage et al., 2012; Carvajal et al., 2017; Maier et al., 2020). Described as either “train of scours” or “erosional steps” (e.g., Fildani et al., 2006, 2013), these features have been related to the upstream migration of cyclic steps sculpted by turbidity currents, which in essence are a train of spaced hydraulic jumps that in modern systems have been shown to occur where sediment gravity flows lose confinement, like in the overbank area of an adjacent active channel (e.g., Fildani et al., 2006; Armitage et al., 2012) and downflow of a break in slope (Carvajal et al., 2017; Maier et al., 2020). Here, regardless of the process that causes the avulsion or where it occurs, it is postulated that erosional knickpoint retreat not only charged the flows with fine-grained sediment, but in conjunction with the depositional bulge created by the underlying stratal unit, drove the systematic back-stepping of successive stratal units (e.g., Fig. 14B) (see also Carvajal et al., 2017). Moreover, in the case of the CF example, the knickpoint scour is not part of a train of scours but instead a solitary feature at the terminus of the channel, and like in the WSG, depositional conditions were controlled by changes in effective fluid viscosity caused by the incorporation of eroded mud rather than changes in hydraulic flow conditions and/or transformations in flow type (Ningthoujam et al., 2022).

A similar history of knickpoint formation and headward migration associated with local submarine channel avulsion has been interpreted from seismic images in the Late Pleistocene Amazon and Rhône submarine fans by Pirmez et al. (2000). Downflow of the avulsion node, these authors noted a lobate, high-amplitude (i.e., sand-rich), sheetlike deposit they termed a “high-amplitude reflection packet” (HARP), which here, like in Terlaky and Arnott (2014), is termed an avulsion splay. As stratal units accumulate, local topography increases and eventually causes flows to be diverted into adjacent topographic lows, such as the thin, fine-grained fringes of the depositional continuum. With these diverted flows experiencing an abrupt gradient change, a new knickpoint is created (see also Pirmez et al., 1997), and the process repeats with the deposition of the next stratal unit of matrix-rich lithologies. This combination of knickpoint migration and local topographical steering continues to form matrix-rich deposits until the supply of mud from local seafloor erosion becomes exhausted, quite possibly due to scour down to a sufficiently compacted subsurface layer, the avulsion node healing, the muddy substrate becoming buried in sand, or the parent channel diverting to a new location.

Matrix-rich and associated matrix-poor sandstones typically underlie sand-rich stratal elements like slope channel, distributary channel, and terminal splay deposits in the Windermere Supergroup (WSG) and terminal splays in the Cloridorme Formation (CF). This association is interpreted to be related to the activation of the local transport system caused by an upflow avulsion. Avulsion formed a wall jet that locally scoured the mud-rich inter-channel area and charged the flow with fine-grained sediment, namely very fine sand to clay, in addition to low-density mud clasts. The local incorporation of easily suspended fine-grained sediment resulted in an abrupt increase in local effective viscosity and a dramatic change in the characteristics and reduction in the diffusion of turbulence required to maintain the particle suspension. This caused the sediment suspension to rapidly exceed transport capacity and collapse, forming a negligibly sheared suspension that deposited a systematic along-flow depositional continuum consisting of matrix-poor sandstone (MPS) to muddy sandstone (MS) to bipartite facies (BF) and then sandy mudstone (SM) over a distance of hundreds of meters and collectively overlain by a thin-bedded, traction-structured sandstone and/or silty mudstone cap. Deposition of this transect occurred along the margins of the high-energy wall jet in the WSG, but at the downflow terminus of the flow in the CF.

Vertically, and at the bed scale, similar lithofacies tend to preferentially stack and build up multi-bed stratal units two to nine beds thick. The grain-size distribution within a multi-bed stratal unit shows a subtle but nevertheless consistent upward-fining trend, whereas the matrix content changes little upward. Additionally, within a discrete multi-bed stratal unit, individual beds undergo a transformation from one facies to the next in the lateral depositional continuum at about the same lateral position. Collectively, these observations suggest that beds within a discrete multi-bed stratal unit are more or less texturally self-similar and therefore deposited from flows with similar hydraulic and compositional character, albeit exhibiting a minor but nevertheless progressive upward decrease in flow energy. This spatial and temporal regularity is argued to have been generated by multiple avulsion-related surges of a single, gradually waning turbidity current. Importantly, the time between the successive surges is estimated to have been of sufficient duration to allow for the emplacement of the traction-structured sandstone and overlying mud cap. Furthermore, the upward stacking at successive stratal units is dominated by side-stepping with minor back-stepping in the WSG compared with mostly back-stepping in the CF. Overall, the initiation, back-stepping, and side-stepping of stratal units are interpreted to be related to a combination of knickpoint migration and local topographical steering of the flows, which continued until the supply of mud from local seafloor erosion became exhausted.

We thank journal reviewer Marco Patacci, an anonymous reviewer, and Associate Editor Andrea Fildani for their helpful suggestions and insights on an earlier version of this manuscript. Funding for this research has been generously provided by members of the Windermere Consortium (Occidental Petroleum [OXY] and Cenovus [Husky Energy]), Husky Energy Fellowships in Sedimentology and Petroleum Geology to Ningthoujam and Wearmouth, and a Natural Sciences and Engineering Research Council of Canada (NSERC) Discovery Grant to Arnott.