Large-scale modification of submarine geomorphic features on the Cascadia accretionary wedge caused by catastrophic flooding events Open Access

Submarine canyon formation along continental margins results from the interplay of retrogressive erosional processes, deposition, and levee formation (Orange, 1999; Lamb et al., 2008; Willems et al., 2011; Bourget et al., 2011). Typically, erosion in terrestrial systems is translated to deposition in the deep sea through erosion of the continent and sediment transport through river systems to coast and shelf environments and down submarine canyons and channels to the deep sea.

Active margins have tectonic forces that need to be considered when investigating marine depositional systems (Orange, 1999; Lamb et al., 2008; Mountjoy et al., 2009; Mouchot et al., 2010). Accretion of sediments from subduction creates fold-and-thrust belts that, along with sediment transport from fluvial systems, are subject to coeval submarine canyon incision. Tectonically over-steepened slopes increase landslide occurrence, which, with repeated large seismic events, dislodge and transport large volumes of sediment to the deep sea (Adams, 1990; McAdoo et al., 2000). In rare cases, erosional processes translate across terrestrial-marine boundaries, eroding terrestrial systems as well as deep-sea environments (Gupta et al., 2007; Willems et al., 2011). It is clear that large high-energy events can have profound impacts on regional geomorphology (Baker, 1973; Clifton, 1988; Poag, 1997; Baker and Kale, 1998; Dawson and Stewart, 2007). These convulsive events can provide the energy needed to significantly shape, reform, and alter sediment routing paths within a sedimentary system (Gupta et al., 2007).

In this study, we use recently collected high-resolution data sets offshore of Washington State, USA, to describe and interpret features on the seafloor¹. Through our careful observations and interpretations, we build a hypothesis for what event(s) may have contributed to the geomorphology we image today.

The Cascadia margin is characterized by the subduction of the Juan de Fuca plate beneath the North American plate at a convergence rate of 34.2 ± 1.2 mm/yr near the latitude of Seattle, Washington (DeMets et al., 2010). This convergence coupled with high sediment supply has created a wide fold-and-thrust belt above the subduction zone (Davis and Hyndman, 1989). A 3–4-km-thick Pleistocene sediment package deposited on the abyssal plain has been folded and thrust into a 30–50-km-wide accretionary wedge offshore Washington and Oregon (Davis and Hyndman, 1989; Goldfinger et al., 1992, 2003; Westbrook, 1994).

The Cascadia margin contains several large submarine canyon systems that connected to fluvial systems during Pleistocene sea-level lowstands (Nelson, 1976). This connection created pathways for sediment to bypass the shelf and be transported to the deep sea (Nelson, 1976; Goldfinger et al., 2012; Normark and Reid, 2003) (Fig. 1B). In contrast, Pleistocene deposition during sea-level highstands occurred mainly in shelf basins and upper canyon heads (Posamentier et al., 1988; Goldfinger et al., 1992; McNeill et al., 2000).

Offshore of Washington, Nitinat Fan, originally described by Carson et al. (1974) and Normark (1970), is fed by Nitinat and Juan de Fuca Canyons, which coalesce near the base of the slope (Fig. 1B). Below the confluence, the Juan de Fuca channel turns south, cutting off older channels on the middle to outer fan (Carson et al., 1974; Barnard, 1978). The Juan de Fuca channel then parallels the base of the slope and turns southwest to merge with Willapa Deep-Sea Channel (Carson et al., 1974). During the present sea-level highstand the Pleistocene canyon systems are largely relict except for ∼400–500 yr recurrences of seismically induced flushing of canyon systems (Carpenter et al., 1982; Adams, 1990; Goldfinger et al., 2012).

Astoria Canyon is the largest submarine canyon offshore of Oregon, its canyon head ∼20 km west of the mouth of the Columbia River. The width of Astoria Canyon (rim to rim) is between 2 and 11.5 km; the width of the canyon floor is ∼2 km and ranges from 0.5 to 5 km (Carlson, 1968). In the middle section of the Astoria Fan, Astoria Canyon splits into two channels. One channel continues south paralleling the deformation front, the other, farther offshore, to the southwest (Carlson, 1968; Nelson, 1976; Goldfinger et al., 2003, 2012).

Willapa Canyon, the main canyon for this study, lies north of Astoria Canyon and marks the boundary between Nitinat and Astoria submarine fans (Fig. 1C). Grays, Quinault, Guide, and Willapa Canyons merge at ∼1600 m water depths to form the larger Willapa Canyon (Fig. 1B). Sections of Willapa Canyon have a wide range of channel widths ranging from >1 km to <7 km (Royse, 1964; Carlson, 1968; Goldfinger et al., 2003, 2012). The Willapa system lacks a major depositional fan, which suggests lesser sedimentation there compared to the larger Astoria and Juan de Fuca systems.

We created high-resolution (15 × 15 m) bathymetric grids covering part of the Cascadia accretionary wedge (Fig. 1C) using multibeam bathymetric data collected for the Cascadia Open-Access Seismic Transects (COAST) and Cascadia Initiative projects (Holbrook et al., 2012; Toomey et al., 2014). Multibeam bathymetric and backscatter data were collected on the R/V Thomas G. Thompson using a Kongsberg EM302 multibeam system, and on the R/V Marcus G. Langseth using a Kongsberg EM122 multibeam system. Both ships used high-precision differential global positioning system and real-time attitude corrections. Expendable bathymetric thermographs were deployed twice daily, at minimum, to record and correct for water-column sound-velocity variations.

Geological Long-Range Inclined Asdic (GLORIA) backscatter data collected by the U.S. Geological Survey during the 1980s and 1990s were also used for seafloor and near-seafloor composition interpretations (EEZ SCAN 84 Scientific Staff, 1986). Our experience suggests that penetration can range from a few meters in sandy substrates to up to 30 m in muddy or silty environments, thus these data are imaging an amalgamated time period rather than simply imaging the most recent surface materials.

A hull-mounted Knudsen 3260 sub-bottom profiler, with a frequency sweep of 2–6 kHz, recorded three seismic profiles crossing two frontal thrust erosional features on the R/V Thomas T. Thompson cruise TN265 in 2011. Processing of these data included application of a bandpass filter and a time-varying gain. Multichannel seismic (MCS) data were acquired on the R/V Marcus Langseth with 36 airguns in a 6600 cu. in array which provided an acoustic signal recorded on a 636-channel, 8 km streamer sampling the subsurface every 6.25 m (Holbrook et al., 2012). These data were pre-stack depth migrated following conventional seismic processing steps including: bandpass filtering, trace editing, source deconvolution, multiple elimination, velocity analysis, and horizon-based tomographic velocity updates (Yilmaz, 2001).

The high-resolution data presented here allow a new analysis and interpretation of geomorphic features in the Cascadia accretionary wedge. The three primary features (A, B, and C) studied are unique geomorphic features where the unknown underlying mechanism prompted this investigation. These features consist of (1) steep (between 18° to 25°) semi-circular headwall scarps, (2) channel-like troughs that extend from the headwall scarps through the frontal anticlinal ridge, (3) margin-parallel troughs in the most landward piggyback basin that extend back from the headwall scarps to the Willapa Deep-Sea Channel, and (4) large depressions, west to west-northwest trending, that are partially sediment filled and extend west from eroded ridge sections by as much as 16 km. Previous investigations of sediment routing, using coarser bathymetric grids, interpreted these features as inactive channels in the Cascadia basin-channel systems (Goldfinger et al., 2012).

Geomorphic feature A is the farthest north, located where the frontal thrust ridge abruptly steps to the west by ∼3.6 km (Fig. 2). The feature incises the frontal thrust ridge and the flanking margin-parallel structural trough (piggyback basin) to the east, terminating in a steep (22°), 165-m-high headwall scarp. This trough appears to be a channel bounded by subtle 15-m-high parallel walls that connect to feature B to the south. This channel is being filled with sediment as indicated by very low backscatter response relative to the bounding levees and the crosscutting Juan de Fuca Canyon levee. A few scattered landslide debris blocks are visible and derived from the anticlinal ridge on the channel’s eastern flank (Fig 2).

The headwall scarp bounds a semi-circular depression ∼3.5–4 km in diameter within the structural trough, with a 1.5–2-km-wide neck that extends west through the frontal thrust ridge (Fig. 2). The multibeam bathymetry data image several small landslide blocks within the semi-circular depression, notably a large (0.3 km²) block that sits ∼90 m above the depression floor. The multibeam bathymetry and MCS line each image a ∼60-m-high west-side-up fault scarp on the west flank of the depression (Figs. 2, 3B). The scarp can be traced north and south beyond the sidewalls of feature A, as confirmed by seismic profiles (Fig. 3B), and is interpreted as the landward-vergent scarp of the frontal thrust fault of the Cascadia accretionary prism.

A 6.4-km-wide circular depression in the abyssal plain occurs at the base of the excavated frontal ridge. The deepest part of this depression is near circular and ∼43 m lower than the surrounding abyssal plain. Seismic-reflection data (Figs. 3B and 3C) reveal that the depression is underlain by a channel filled by as much as 160 m of sediment (assuming a seismic velocity of 1600 m/s). The western margin of this depression is 6 km from the frontal ridge and ∼7 km east of the active Juan de Fuca channel (Goldfinger et al., 2012). The Juan de Fuca channel and its most recent levees show no discernable disruption in the new bathymetry. A chirp sub-bottom seismic profile images the upper ∼50 m of the circular depression on the abyssal plain (Fig. 3A). This profile shows flat-lying low-amplitude reflections that onlap the walls of the scour depression and thicken toward its center.

About 2 km north of feature A, the new bathymetry images a second excavated section of the frontal anticline (Fig. 2). This excavated section is between 0.5 and 1.5 km wide and isolates a piece of anticline ridge between this embayed section and the main one at feature A.

Geomorphic feature B is lies ∼18 km south of feature A where the frontal thrust ridge changes strike from 342° to 326° (Fig. 2). It incises 1.7 km east of the north-south frontal anticline crest into the frontal piggyback basin terminating in a steep (25°) 278-m-high scarp (Fig. 2; Fig. 4, profile B-B′). The headwall scarp connects the excavated section to a 2-km-wide, 15-km-long trough. This trough is apparently a channel that crosses between the frontal and second piggyback basins from feature A southward toward the modern Willapa Deep-Sea Channel (Fig. 2). At the connection between the trough or channel and Willapa Canyon, a ∼125-m-high (crest to thalweg) barrier paralleling Willapa Canyon separates the excavated ridge section from the active Willapa system (Fig. 2; Fig. 4, profile B-B′). This height includes the step-up from the Willapa thalweg into the piggyback basin (∼90 m) which is topped by a ∼35-m-high levee paralleling the main Willapa Canyon (Fig. 4, profile B-B′). This levee height above the trough floor is similar to the downstream average levee heights of ∼27 m on the abyssal plain described by Griggs and Kulm (1970) and verified by modern data (this study).

Feature B has cut back into the adjacent piggyback basin with a 1.3–2-km-wide stair-stepping channel to the abyssal plain. The channel extends away from the excavated ridge section and reaches depths ∼100 m below the adjacent abyssal plain. The feature B depression is more elongated than that of feature A, extending west at least twice its width to an intersection with the Holocene Juan de Fuca channel (Fig. 2).

This elongate depression is imaged in both the bathymetry and on north-trending sub-bottom profiles (Figs. 5A and 5B). This likely paleochannel is ∼2 km wide at the mouth of feature B and widens to 4.6 km at the location of Figure 5A (see Fig. 1) where ∼17 m of sediment fill the channel. The paleochannel trends 270° which is at nearly right angles to the present local slope and the modern Juan de Fuca channel.

Geomorphic feature C is the smallest (2 km wide by 3.5 km long) of the three excavated sections on the frontal ridge. Its headwall scarp is less steep (18°) than those of the other features, with its excavated section not extending completely across the frontal anticline into the piggyback basin (Fig. 2).

A depression as deep as 29 m at the base of the excavated frontal ridge extends west ∼3.5 km across the surrounding abyssal plain (Fig. 2). Four kilometers south of feature C, a secondary paleochannel leads around the flank of the frontal anticline along a probable fault in the wedge (Fig. 2).

The toe of the accretionary wedge has numerous submarine landslides of varying sizes (Fig. 6). Two-dimensional depth profiles E-E′ and F-F′ (Fig. 4) cross two landslide scars on the frontal ridge that do not extend past the frontal ridge peak. Most of these slides are tabular bedding-plane slides from the frontal limb of the landward-vergent ridge, which slopes west an average 6° to 7°. Most of the slide debris appears to be tabular blocks that disaggregated from larger slide blocks during landslides (Fig. 6).

Debris piles of a few submarine landslides have runout of as much as ∼5–6 km from the base of the landslide scarp (Fig. 6). The landside debris deposits are typically partly buried by abyssal-plain sediments but include geomorphically “fresh” landslide blocks sitting as much as 85 m above the abyssal plain. None of the imaged landslide debris deposits are associated with large excavated ridge sections or elongate scour depressions.

The upper reach of Willapa Canyon has an active channel system (400 to 500 m wide), bounded by small levees, that is much smaller in comparison to the much wider (7 km) host channel (Fig. 7, profiles A and B). Profiles A and B on Figure 7 also show paleochannel walls ∼300 m above the active channel floor.

The lower reach of Willapa Deep-Sea Channel, as imaged, crosses and has eroded through several anticlinal ridges (Fig. 7). Profiles C and D in Figure 7 highlight the wide Willapa paleodrainage and its narrower inset Holocene channel. At the lower reach, the paleodrainage width shrinks to 2.3 to 3.5 km while the inset active channel widens to ∼1 km. It is also at this section of canyon that the smaller Quinault Canyon system merges, between profiles C and D in Figure 7. These bathymetric profiles show a secondary drainage north of the main Willapa drainage (Fig. 7, profiles C and D) characterized by eroded anticlines and depressed bathymetry, starting west of profile B in Figure 7.

The wide (4 to 8 km) paleodrainage system associated with the active Willapa submarine channel suggests an underfit system. Underfit systems are defined as drainage areas too small to correlate with current channel characteristics. While it would not be unusual for the Holocene channel to be underfit relative to the Pleistocene channel, the size difference in this case is strikingly exaggerated in comparison to other Cascadia systems. For example, the Pleistocene Astoria system was connected to the Columbia River and is the primary drainage system for western North America. Astoria Canyon’s Pleistocene canyon is 2–3 km in width with the main Holocene channel ∼500–700 m in width where observable. In contrast, the width of the Pleistocene Willapa system at ∼5 km is more than seven times that of its ∼700 m Holocene inset channel. This large paleodrainage could potentially be explained by increased flow down Willapa Deep-Sea Channel during the Pleistocene as is the case for other systems throughout the Cascadia margin. However, it was primarily the Astoria system that received the drainage from the Columbia River during the Pleistocene, leaving Willapa Canyon with no discernable connections to other major river systems, only the much smaller adjacent coastal drainages.

The ages of the three unique geomorphic features identified in this study (features A, B, and C) can be indirectly constrained in three ways: (1) overlap and crosscutting relationships with the Juan de Fuca channel, Willapa Deep-Sea Channel, and scour depressions, (2) the size and location of post-erosion event(s) levees, and (3) the location and offset of the features by the frontal thrust fault identified at feature A.

The multibeam bathymetry and sub-bottom profiles show that the scour associated with feature B extends from the toe of the deformation front out to, and likely past, the active Juan de Fuca channel that presently overlaps it (Figs. 1C and 2). The extent of feature B’s scour into the abyssal plain, both in its basal depth of 100 m below the adjacent plain near the deformation front and its 17 km minimum extent to the west, indicates that the Juan de Fuca channel would have received the erosional currents that formed feature B if the two systems were coeval. But now the Juan de Fuca channel path adjacent to feature B is smooth and undisturbed, indicating that the modern channel and levee system developed and evolved after the apparent erosional event(s) that formed feature B. Karl et al. (1989), using GLORIA sidescan data, suggested that the Juan de Fuca channel is young because it can be traced to the canyon directly, unlike other channels on the Nitinat Fan. Griggs and Kulm (1970) came to a similar conclusion because seismic-reflection profiles show little development of the Juan de Fuca channel in the subsurface. From radiocarbon dating of turbidites, Goldfinger et al. (2012) demonstrated that the Juan de Fuca channel and its levees are mostly Holocene features.

Older paleochannels associated with a meandering Juan de Fuca channel show in the high-resolution bathymetry west of the present Juan de Fuca channel except at the latitude of feature B where they are not discernable (Figs. 1C and 2). The area where these paleochannels are absent is on trend with the feature B scour channel, and we suggest that the currents generated from feature B also eroded away these channels. It would have taken at least ∼20 m of erosion to erase the missing section of paleochannel. The Juan de Fuca channel, the main sediment routing path at present, would have had to reform after the feature B erosive event(s).

Feature A’s abyssal plain scour depression neither erodes nor impinges on the active Juan de Fuca channel, nor does it appear to affect the paleochannels farther west. Nothing in the bathymetry or backscatter data suggest a barrier to such erosive effects were they present during the time of formation of feature A. The scour depression of feature A is different from that of feature B; it is more circular and less asymmetrical. The lack of erosion as far west as in feature B, and the more symmetrical depression morphology, suggest that feature A may have formed with less of a horizontal component of erosion during its formation. These observations lead to two possible scenarios: (1) event(s) that created feature A are older than both the most recent activity and levee construction of the Juan de Fuca channel and the western paleochannels, or (2) feature A’s erosion was less robust and was limited to the area near the deformation front, not reaching past the present Juan de Fuca channel to the area of the western paleochannels.

In GLORIA data, feature B’s scour channel is imaged as a 2.3-km-wide zone of dark reflectance that extends from the deformation front out to the active Juan de Fuca channel and likely beyond (Fig. 8). The dark reflectance indicates either a softer substrate or one lacking in specular reflecting particles relative to nearby high-reflectance areas. This suggests a muddier average substrate, mostly likely with smaller particles within the scour and near the surface. The flanks of this scour channel have brighter reflectance, which we interpret as paleolevees built while this scour channel was active (Fig. 8).

GLORIA data also reveal the paleochannel connecting feature B with the active Willapa Deep-Sea Channel via the first piggyback basin (Fig. 8). This channel is also characterized by low reflectance and corresponds with the trough imaged in the high-resolution multibeam bathymetry. We interpret this low-reflectance zone as fine-grained sediments similar to those of the scour channels, indicative of a paleochannel or low-activity channel.

A ∼20-m-high levee crosses the axis of the paleochannel on the north flank of feature B in the piggyback basin east of the frontal ridge (Fig. 4, profile A-A′). This levee suggests that feature A was created before or coeval with feature B, then a levee was built when most of the flow occurred in feature B (Fig. 4, profile B-B′). Another levee in this paleochannel presently impedes the connection between features A, B, and C and the active Willapa Deep-Sea Channel. This suggests that the next progression was a shut-off of feature B with a shift of most or all flow down the now-active Willapa Deep-Sea Channel. This levee sits ∼100 m above the thalweg of Willapa Deep-Sea Channel and is ∼30–35 m above the paleochannel floor. Considering that this levee crosscuts the entrance to the paleochannel and parallels the main Willapa Canyon, the levee formed near or after the time when flow into features A, B, and C ended and most or all flow was directed down Willapa Deep-Sea Channel. Given that Holocene fill rates of the nearby Juan de Fuca and Willapa channels and levees in the Holocene range from 30 to 60 cm/1000 yr, (Goldfinger et al., 2012, 2017), building the levee would have required 50–112 k.y. This indicates that either rates of levee construction would have had to have been much larger than Holocene rates of sedimentation, or the levees are much older than late Pleistocene.

The relationships indicated by the crosscutting of feature B by the Juan de Fuca channel and its levees would place a relative age of feature B as older than the eroded paleochannels, but younger than the latest Juan de Fuca channel incision. The ages of the older paleochannels are presently unconstrained.

The ∼60-m-high fault scarp in the frontal ridge at feature A is inferred to have resulted from slip on the frontal thrust fault in the accretionary wedge after the erosion event that created feature A. There is no apparent incipient thrust fault imaged in the MCS profile that could be accommodating accretionary wedge deformation farther west (Fig. 3B). The scarp height and estimated fault dip (38°) from the MCS section suggests ∼76 m of fault slip, a relatively small amount compared to the ∼300 to 400 m of Holocene shortening expected along the subduction zone based on plate convergence rates (DeMets et al., 2010). This relatively small estimate for fault slip after feature A creation suggests (1) a very young age for the fault itself, (2) a young age of the offsets seen in features A and B, and therefore a young fault age, or (3) a slow slip rate, less than plate convergence at the frontal thrust.

A good candidate for flows far larger than expected in the Pleistocene are the massive glacial jökulhlaups known as the Missoula flood events (Bretz, 1928; Baker, 1973). Landforms across eastern Washington attest to the massive erosive power of the Missoula events which have shaped the primary Quaternary geomorphology of eastern Washington and the Columbia River Gorge (Bretz, 1928; Patton and Baker, 1978; Benito and O’Connor, 2003; Peterson et al., 2011). One of the most visually striking examples is the numerous elongate ravines with steep headwall scarps eroded into Columbia River Basalt, known as “coulees” (Baker, 1973). Based on similarity to terrestrial geomorphic features created by the same mechanism, we see it fitting to consider the submarine erosional features described in this study as the first observation and recognition of “submarine coulees”.

We further suggest that these events severely eroded several anticlines which now terminate abruptly at the boundaries of the broader drainage as surficial features. These anticlines, along with a broad section of the accretionary wedge around Willapa Deep-Sea Channel, have the characteristic of a broadly eroded zone, with no other obvious explanation for the truncated folds. The piggyback basins to the east of features A and B both show bathymetry significantly lower than other piggyback basins to the north. We suggest that this entire section experienced extensive erosion for a short period of time.

Flow rates of the largest Missoula floods have been estimated to be over a thousand times that of the modern (∼6970 m³/s) average Columbia River discharge (Naik and Jay, 2005). The Missoula floods primarily emptied through the Columbia River and presumably down the Astoria and Willapa submarine channels out to the deep sea where Missoula flood deposits have been identified >1000 km away in Escanaba Trough (Brunner et al., 1999; Zuffa et al., 2000; Normark and Reid, 2003). During sea-level lowstands, fluvial systems connected to submarine canyon systems which allow sediment to move directly to the deep sea. The more tightly coupled systems likely allowed the deep-sea environment to be more susceptible to alteration by large flooding events.

The Missoula floods occurred between ca. 19,500 and 15,000 yr B.P. (Brunner et al., 1999; Zuffa et al., 2000; Benito and O’Connor, 2003; Waitt, 2016). Normark and Reid (2003) and Zuffa et al. (2000) suggested that the Missoula flood outpouring from the Columbia formed large hyperpycnally generated turbidity flows down both Astoria and Willapa channels. The inferred flow paths and the limited timing constraints on the origin of geomorphic features A, B, and C are in broad agreement with this, and we therefore conclude that the Missoula floods are a likely mechanism for the erosion seen at this section of the accretionary wedge.

We suggest that turbidity flows associated with the Missoula floods flowed down Willapa Canyon and encountered the back side of the first frontal ridge that deflects Willapa Canyon south by nearly 90°. The abrupt blockage ∼500 m high turned these flows south in the main channel and north into the hanging valley of the piggyback basin (Fig. 8). This would have been the second-least-resistive path and a logical overflow path for such events that could not overtop the frontal anticline. We propose that the high-volume and high-density currents made this jump into the piggyback basins to the north and to the south and subsequently eroded through the frontal anticline, concentrated at weak and/or low points along the inboard part of the frontal anticline.

The Missoula floods were not a single catastrophic event that occurred near the end of the Pleistocene (Atwater, 1986; Waitt, 1985, 2016; Benito and O’Connor, 2003). Between 90 and 100 separate flooding events of various sizes occurred on average every 40 yr starting ca. 19,500 yr B.P. (Benito and O’Connor, 2003). The early flooding events are thought to have been the largest and most catastrophic, with the size and intensity of flooding waning near the end of the Pleistocene (Waitt, 1985, 2016; Benito and O’Connor, 2003). The diminishing of flow size may explain the thickness of sediment fill imaged at feature A’s plunge pool (∼165 m). The diminishing flows over time would have filled earlier excavations from the largest flows. The chirp profile in Figure 3 indeed shows decreasing bed thickness toward the top of the section. The recognition of these later, smaller-sized Missoula floods also offers an explanation for the 35-m-high levee blocking the paleochannel and Willapa Canyon. Smaller floods would have still sent large amounts of sediment down Willapa Deep-Sea Channel. These smaller floods could have been more depositional than erosive events, and levees for the smaller events could have retreated to the top of the main Willapa Canyon sidewall as observed in the present bathymetry.

With our interpretation that the geomorphic features identified in this study were largely formed during abnormally large sediment flows, we explore how the hydrodynamics and sediment dispersal of large-volume turbidity currents would evolve down the Willapa drainage system if the bathymetry at features A, B, and C were restored. We used a three-dimensional turbidity-current modeling module in Move 2015 software package (https://www.mve.com/), which is based on turbidity-flow physics (Waltham et al., 2008). Such modeling, used for hydrocarbon exploration of turbidite-filled basins, is easily adaptable for qualitatively exploring how large-volume flows evolve in this paleodrainage system.

Waltham et al. (2008) reported that turbidity-current modeling can be approached using forward modeling. In a forward model, flow-distribution patterns can be predicted from input of flow parameters and particle information. In the case of modeling a conceptual Missoula flood–style flow, where we lack any physical deposit samples, we used reasonable estimates of flow characteristics to explore how increasing flow volumes and concentrations would move through the paleobathymetry.

We reconstructed the bathymetry that likely existed before the erosive event(s) created features A, B, and C. This reconstruction was done by extending the height and geometry of the frontal ridge across the erosion zones (Fig. 9). We used flow parameters that were consistent with literature on physical properties of turbidity flows (Stow and Bowen, 1980; Mulder and Alexander, 2001; Manville and White, 2003; Egawa et al., 2013). We used sediment concentration values ranging from 1% to the upper limit of 9% for suspended flows (Bagnold, 1962; Mulder and Alexander, 2001). Inflow heights ranged from 50 to 300 m, and inflow width was estimated by the canyon width of the inflow regions.

Modeling indicated that a minimum volume flow of >15 km³ of fine to medium sand with concentrations of 3.5% initiated at the head of Willapa Canyon has the potential to fill and inundate the entire paleodrainage system along with the piggyback basin (Fig. 10D). Flow initially fills the upper regions of the canyon system, followed by filling the paleodrainage and secondary systems (Fig. 10C). Continued input into the system inundates the main drainage system along with entering and filling of the piggyback basin between the first two thrust ridges (Fig. 10D). With the reconstructed bathymetry lacking a flow path out of the basin, flow waters ultimately filling the basin and overtop a low point on the ridge. The low spot in this case is a landslide scar where a section of ridge top has been removed, allowing for basin flow to drain. While this model run only simulates a potential large-flow scenario, it shows the plausibility of our conceptual model of formation of features A, B, and C.

Additionally, we considered the seismic profile within feature A’s depression, which shows seismic units between 1 and 3 m thick (Fig. 3A). The flow volumes needed to deposit beds of that thickness within feature A’s depression are at least 17 km³. This model result further favors the need for large flows, potentially the later Missoula floods, to have entered the Willapa system, accounting for deposition within feature A’s depression.

Discharge rates of the largest Missoula floods could easily have produced inflow volumes >15 km³. Using the largest discharge rate of 6 × 10⁶ m³/s over the course of a single day would produce flow volumes >500 km³. Estimating a flow concentration of 3.5% would equate to ∼18 km³ of inflow deposits (a parameter used for modeling). This calculation for one day of this flow is likely a minimum volume, with the estimates from previous studies suggesting that Missoula flooding could have taken multiple days to drain (Denlinger and O’Connell, 2010).

We present new high-resolution multibeam bathymetry and seismic-reflection data collected offshore Washington on the Cascadia accretionary wedge and abyssal plain. These data more clearly than previous datasets reveal submarine erosional features on the toe of the deformation front, the accretionary wedge, and the abyssal plain.

Three features (A, B, and C) connect to a paleochannel that was sourced from overbanking of the main Willapa Canyon. This confluence valley is currently overlapped by a 35-m-high levee paralleling Willapa Canyon. This levee is >130 m above the thalweg of Willapa Canyon and was likely formed subsequent to the occupation of the paleochannel by late Pleistocene–Holocene turbidity currents.

These features represent three areas of the backlimb of the frontal thrust anticline that have been deeply incised and excavated back into the adjacent piggyback basin. The morphology of these features is a poor match to typical tabular blocky slides, common on the Cascadia margin. They lack blocky debris, or long vertical drops and wave fields typical of landslides and plunge pools. Rather, they extend well into the landward piggyback basin, and two of them have asymmetrical scours extending deep into the abyssal plain and across the regional slopes. These features are a better morphological fit to large erosional features associated with the Missoula flooding events, features known as “coulees,” and should be regarded as “submarine coulees” based on their morphology and mechanism of formation.

Model results suggest the plausibility that the erosional features (A, B, and C), abyssal plain, and paleodrainage system of the Willapa submarine canyon were formed by the Pleistocene Missoula floods. Input flow volumes to occupy, much less create and erode, the paleochannel, features A, B, and C, and oblong abyssal plain scour features are at a minimum about one order of magnitude greater than Holocene flows that can be fit to core constraints in the Willapa Canyon drainage. These results give insights into the influence that glacial outburst floods and perhaps other catastrophic phenomena can have on the marine environment evolving more slowly through other processes.

The potential for a large section of the accretionary wedge to be profoundly altered by catastrophic outburst floods highlights the need to consider these events when doing geomorphic interpretations in other mid- to high-latitude submarine environments. While the framework and general morphology of an active margin will be dictated by the long-lived interplay of tectonic forces and sediment transport and deposition, convulsive events can and do play a role in shaping continental margin architecture.

New multibeam and seismic-reflection data were collected under National Science Foundation grants OCE-1150628 and OCE-1137986. Thanks also go to the crew and science party of the R/V Langseth and R/V Thompson during the research cruises. Many thanks go to Sam Johnson for his detailed reviews and edits. We would like to acknowledge and thank R.B. Wyatt and an anonymous reviewer. Their thoughtful comments and suggestions greatly improved the manuscript. Paradigm’s Focus and GeoDepth software packages were used for multi-channel seismic-reflection data processing. Midland Valley’s sediment modeling package (Move) was used for turbidity-current simulations.