Deformed Latest Pleistocene Fluvial Terraces Reveal Complex Active Faulting within Tectonic Transition Zone, Mendocino Triple Junction, Northern California

The Mendocino Triple Junction (MTJ) has been the archetypal tectonic setting to understand overprinting of various tectonic processes, triple junction behavior (specifically transform–transform–trench type triple junctions), and the signal of plate coupling in subduction zones (e.g., Furlong and Govers, 1999; Furlong and Schwartz, 2004; McKenzie and Furlong, 2021). The MTJ is defined by the region where the San Andreas fault zone (SAF), the Mendocino fracture zone, and southern extent of the Cascadia subduction zone meet (Fig. 1a) (Kelsey and Carver, 1988). This complex zone experiences northeast‐directed contraction of the southern CSZ, transitioning to north‐northwest‐directed translation within the broad SAF transform margin to the south. The MTJ migrates north‐northwest through time, proceeded by crustal shortening and accompanying crustal thickening and uplift driven by Mendocino crustal conveyor (MCC; Furlong and Govers, 1999; Lock et al., 2006; McKenzie and Furlong, 2021).

Presence of SAF‐type deformation north‐northeast of the northernmost mapped extent of the SAF zone has been suggested through geologic and geodetic data (Kelsey and Carver, 1988; Williams et al., 2006; McKenzie and Furlong, 2021). Because the geological and plate‐velocity representation of these two tectonic regimes exhibit different spatial extents, the transition between them remains enigmatic (Kelsey and Carver, 1988; Williams et al., 2006). However, structural accommodation of northward translational deformation associated with the SAF has been recognized to include east–west‐trending folds in the Eel River Valley, oblique slip on the Little Salmon fault (LSF), Mad River fault zone (MRFz), and Russ fault (Kelsey and Carver, 1988; Williams et al., 2006), and a possible component of dextral motion along the Goose Lake fault (Ladinsky et al., 2019; McKenzie and Furlong, 2021; Bold, 2022; Fig. 1a). This complex expression of crustal deformation also includes distributed dextral shear strain to the south and southeast via the Maacama fault and the Bartlett Springs/Lake Mountain fault zone (Williams et al., 2006; Leroy et al., 2012; Fig. 1a).

Efforts toward deconvoluting the compound influences on the region surrounding the MTJ have been multidisciplinary and include seismologic (e.g., Dengler and McPherson, 1993; Oppenheimer et al., 1993; Murray et al., 1996; McCrory, 2000), geodetic (e.g., Williams et al., 2006; McKenzie and Furlong, 2021), geomorphologic (e.g., Merritts and Vincent, 1989; Clarke and Carver, 1992; Koehler, 1999; Bennett et al., 2016; Clubb et al., 2020), and paleoseismic approaches (e.g., Carver and Burke, 1988; Hemphill‐Haley and Witter, 2006; Ladinsky et al., 2019). Tectonic activity, expression, and characteristics of minor faults in this area may be key in understanding how crustal deformation manifests in complex tectonic environments with hybrid signals. We build on this work by identifying and characterizing a previously unrecognized fault trace along the southern LSF using tectonic geomorphology and paleoseismic techniques.

In this article, we focus on a newly recognized fault informally named the Van Duzen fault (VDF) located in the Van Duzen River drainage south of the LSF (Fig. 1). The LSF is one of the southernmost active thrust faults within the onshore fold‐and‐thrust belt associated with the CSZ and lies proximal to the transition from compressional to dextral stress across the MTJ where the SAF zone meets the southern extent of the CSZ (Woodward‐Clyde Consultants, 1980; Carver and Burke, 1988; Kelsey and Carver, 1988; Clarke and Carver, 1992; Hemphill‐Haley and Witter, 2006). Thus, this southern zone along the LSF is an ideal location to characterize strain associated with the complex region of transitional stress regimes within the MTJ. This site is particularly advantageous because of (1) the known tectonic features attributed to the southern extent of compression (e.g. LSF); (2) well‐preserved geomorphic markers (young fluvial terraces); and (3) available high‐precision topographic data which reveal geomorphic features otherwise obscured by dense vegetation. The study addresses two main objectives to better characterize the VDF including (1) its geologic slip rate based on an assessment of fluvial terraces and scarp heights as indicators for the relative timing and amount of deformation and (2) the character of deformation through documentation of a paleoseismic trench.

The study area is in the vicinity of the southeastern terminus of the mapped LSF and encompasses the Wolverton Syncline formed in deformed Quaternary/Neogene Carlotta formation sediments (dominantly nonmarine conglomerate, sandstone, and carbonaceous claystone) of the Wildcat Group (Ogle, 1953; Kelsey and Carver, 1988; McLaughlin et al., 2000; Fig. 1). In this area, the strike of the LSF shifts from northwest to west‐northwest and the Carlotta formation forms dip slopes that leave a strong signature on the landscape (Fig. 1). Oswald et al. (2006) interpret right‐oblique transpression accommodation on northwest‐striking faults near the confluence of the Van Duzen River and Root Creek, ∼2 km southwest of the study area. These kinematics are evidenced by fault striae and offsets of fractures and concretion nodules in strath exposures. Lateral shear within this area (e.g., Oswald et al., 2006; Leroy et al., 2012) is thought to represent the northerly migrating dextral shear of the SAF zone in an area dominated by CSZ‐related, northeast‐directed contraction. Therefore, the study area lies within a broad zone of complex upper‐plate deformation—the Mendocino deformation zone (MDZ) of Williams et al. (2006).

Our study focuses on the VDF, a northwest‐trending surface scarp associated with an apparent subvertically dipping fault, which deforms a series of young river terraces south of and subparallel to strike of the main splay of the LSF within the study area (Fig. 1b–d). The fault exhibits a relatively upthrown northeastern block, and traverses several Van Duzen River terrace risers and treads estimated to range from Pleistocene to potentially Holocene in age, as defined subsequently in this study. The fault may represent the eastern extent of the Goose Lake fault, because its linear map expression across river terraces of varying elevation and age suggests similarly steeply dipping geometries (Ladinsky et al., 2019). In addition, the VDF may be associated with the compressional features in the area, the Wolverton syncline and Alton anticline pair, and/or as a potential subsidiary to the LSF.

We characterized flights of fluvial terraces along the Van Duzen River to constrain timing of deformation along the VDF (Fig. 1a,b). Because previous mapping of Van Duzen River terraces has not differentiated detailed flights of river terraces within the study area (Ogle, 1953; Woodward‐Clyde Consultants, 1980), our first step was to do so using available light detection and ranging (lidar) imagery. Each terrace surface along the Van Duzen River was digitized into a polygon.

Next, elevation ranges for each terrace were extracted and made relative to the Van Duzen River thalweg to form groups of appropriate terrace sets, or flights. To achieve these elevation values, we constructed a relative elevation model (REM) by rendering a raster to represent the long profile of the Van Duzen River. This raster was built by extracting elevation values from the Van Duzen River profile, creating polygons that encompass short sections of the river on the short axes and adjacent fluvial terraces on the long axes, applying the river’s intersecting elevation data to each polygon, and converting these polygons with the river elevation data into a raster, the REM. We then subtracted these REM values from those of the original 1 m digital elevation model (DEM) masked to each mapped terrace polygon, which created a normalized DEM of all terrace polygons. From here, we extracted point data from the normalized elevation model and conducted a spatial analysis to assign the proper polygon to suites of point data within each. These point data were then averaged per polygon, leaving each mapped surface polygon with an average relative elevation (including variance) above the Van Duzen River (Table S1, available in the supplemental material to this article).

We used the averaged normalized elevation per individual terrace surface to group terraces into respective flights. We omitted active alluvium (Qal, Fig. 1c) from the dataset, because it was not useful for fault activity analysis. Active alluvium was determined by geomorphic character displayed through remote sensing data as unstable or transient features and plotted equal to or less than 4.3 m above the river channel. The remaining elevation data were plotted in order of ascending topographic value, and the highest/oldest outliers that deviated from the convex trend of the data were removed for subsequent grouping analysis (Fig. S2). The geomorphic reasoning for extracting older terrace data is that with age, the signature of fluvial terraces flights may erode or become ambiguous with the landscape, leaving only a few, or even a single, terrace representative of a flight. These oldest terraces that deviated from the consistent curvature trend of the plotted data were above 30.5 m. The remaining terrace elevation data (>4.3 m and <30.5 m) were grouped into flights using k‐means cluster analysis, in which the number of cluster centroids was defined by the square root of half the count value of input data. We calculated an optimal number of seven centroids and defined each cluster around these centroids (Table S2). Terraces above 30.5 m (relative to the river long profile) were separately split into groups based on a variance envelope derived from the k‐means cluster analysis of the majority data. This variance envelope was defined by the allotted variance per calculated cluster of the >4.3 m and <30.5 m elevation terraces. We defined six more groups of terraces from the >30.5 m data, totaling thirteen groups of terrace flights. Average relative elevation above the Van Duzen River for each terrace group are reported in Table 1.

Finally, we estimate the age of terrace flights, as the age of a tectonically deformed geomorphic or geological feature is central to assess deformation rate. To estimate terrace age, we assume that uplift is equivalent to erosion in the study area, specifically incision of the Van Duzen River. This assumption is validated by measured erosion and uplift rates compared by numerous methods in various areas around the MTJ, which appear to be reasonably consistent (Merritts and Vincent, 1989; Balco et al., 2013; Bennett et al., 2016; Clubb et al., 2020). Regional uplift rates around the MTJ have been established throughout the last few decades to understand the impact of the migrating triple junction, which generally support equating uplift and erosion. For example, Moon et al. (2018) use $B10e$ and $A26l$ to show that catchment averaged erosion rates broadly reflect uplift. Bennett et al. (2016) calculate landslide volume using geomorphic mapping techniques to demonstrate how erosion keeps pace with tectonic uplift through mass wasting processes. To reflect this work and the uncertainties involving the uplift near the VDF, we apply a variety of uplift rates to our model. We use 1.2 mm/yr (Furlong and Govers, 1999; Lock et al., 2006; Bennett et al., 2016), 2.9 mm/yr (Merritts and Bull, 1989; Clubb et al., 2020), and 4 mm/yr (Merritts and Bull, 1989; Merritts and Vincent, 1989). We divide the average terrace height per flight (Table 1) by each of these uplift rates and propagate the error. We then use the mean value of the terrace age estimates using each of these three uplift rates and propagated their error to a single mean value. These mean age values with 1 sigma uncertainty were evaluated using the OxCal modeling program (Bronk Ramsey, 2001, 2008) to perform a Bayesian statistical analysis of age probability distributions (Fig. 2, Fig. S3). Modeled terrace age minima and maxima are reported in Table 1 and indicate a history of terrace formation spanning from ∼43 ka to the present.

We identified a terrace and terrain‐crossing linear scarp using lidar terrain data and selected an accessible area of the scarp along Van Duzen River terrace Qt9 for paleoseismic trench excavation (location shown in Fig. 1b). Here, a shallow, roughly 1.5 m deep, 16 m long exploratory trench was hand‐excavated across the VDF. The trench exposure revealed a series of fine‐grained sediments interpreted to be overbank deposits associated with the modern Van Duzen River atop a stratified, clast‐supported, moderately sorted, well‐rounded, pebble to cobble‐size gravel deposit interpreted to be imbricated fluvial gravels (Fig. 3). This gravel deposit is interpreted to be deformed, because the apparent bedding dips to the southwest into the base of the exposure from the uphill (northeast) portion of the trench excavation, displaying geometry of an antiformal monocline. We were able to trace the top of this gravel deposit through a series of boreholes within the downhill (southwestern) extent of the trench to a depth of ∼2 m demonstrating its continuity from the hanging wall to the footwall (Fig. 3). The architecture of the gravel deposit is best explained by an up‐to‐the‐southwest verging fold within sediments of terrace Qt9 (Fig. 3). We did not observe any faulted material within the trench, likely because it was not excavated deep enough to reach faulted strata, and therefore were not able to construct an event history from this excavation.

The VDF deforms terraces of varying age and elevations with progressive vertical separation on terraces with increased age. The vertical separation measured from high‐precision lidar data along the single identified trace of the fault ranges from 1.3 ± 0.2 m on Qt9, 1.6 ± 0.3 m on Qt8, and 3.7 ± 0.6 m on Qt1 (Fig. 1b,d, Table 1). Vertical slip‐rate estimates along the VDF were determined using measured vertical separation divided by terrace age estimates and range from 0.06 to 0.38 mm/yr (Table 1). These rates do not account for any lateral deformation nor do we estimate dip values for the fault. Quaternary slip‐rate estimates for the VDF are much lower than rates along the LSF, which has an inferred Holocene slip rate of up to 6.1 mm/yr (Carver and Burke, 1988; Carver, 1992). Because fluvial terraces are progressively vertically displaced with age, we assume that they are recording multiple ground‐deforming events.

Though we are able to detail and describe the uncertainty quantitatively in the average vertical slip‐rate estimates, it is more difficult to capture with uncertainty the slip‐per‐event values or a recurrence interval from the data we have presented. Nevertheless, we can make some broad initial approximations. It is reasonable to assume that terraces Qt8 and Qt9 experienced the same number of earthquakes, given their recorded vertical separation values overlap (Table 1). Because terrace Qt1 has approximately double the amount of vertical separation as Qt8 and Qt9, and with the further assumption that one earthquake produced the vertical separation on the younger terraces (Qt8 and Qt9), we may posit that the older terrace (Qt1) is recording two events. With these assumptions, we can roughly estimate 1.2–2.0 m of vertical slip per event. This value does not account for potential lateral offset and includes large uncertainty driven by the lack of stratigraphic evidence per event (i.e., colluvial wedges, offset marker beds, or stratigraphic contacts). It also does not account for potential fault scarp erosion or unrecognized events. Limited to the existing data, terrace age estimates can be applied as rough bounding ages for earthquake events to approach a recurrence interval value. Thus, if two earthquake events are recorded along the VDF scarp, and given the framework of sequential terrace formation, we can reason that the VDF penultimate event occurred after Qt1 formation (but before Qt8 formation) and that the most recent event occurred after Qt9 formation. Though the age estimates for each terrace include very large error, these assumptions suggest a penultimate event within the latest Pleistocene and an Early to Middle Holocene most recent event.

Geodetic evidence for northward migration of the MTJ through transform strain associated with SAF‐style deformation is interpreted to reach north of Humboldt Bay (Kelsey and Carver, 1988; Williams et al., 2006; McKenzie and Furlong, 2021). Quaternary active faults and folds record complex strain patterns that demonstrate the overprinting of competing tectonic forces (McCrory, 2000). Geological evidence for dextral transpression along faults within the study area is noted along the faults at Root Creek mapped by Oswald et al. (2006). The Russ fault and Goose Lake fault expresses evidence for oblique (reverse and dextral) motion to the south and north‐northwest, respectively (Kelsey and Carver, 1988; Ladinsky et al., 2019; Bold, 2022). Further evidence for SAF‐style deformation within the vicinity of the study area is expressed as contractional strain accommodated by relatively east–west‐trending structures (Kelsey and Carver, 1988; McCrory, 2000; Williams et al., 2006).

The west‐northwest orientation of the compressional faults and folds within the study area suggest potential SAF‐parallel shortening at the transition from transform (SAF‐associated) to compressional (CSZ‐associated) tectonic regimes. These structures include the comparatively east–west trend of the LSF at its southeastern extent, the VDF, and the axis of the Wolverton syncline (Fig. 1a). The west‐northwest orientation of these compressional features deviates from the predominant north‐northwest orientation of thrust faults and folds associated with the CSZ zone. We hypothesize that this change in geometry is a manifestation of the leading edge of the northward migration of the MTJ, propagating the SAF system into the CSZ fold‐and‐thrust belt with the rigid Klamath Mountains as a structural back stop (Leroy et al., 2012; McKenzie and Furlong, 2021; Fig. 4). With an orientation evident of SAF‐parallel compression, geological features such as the VDF offer a means to analyze the geological manifestation of transition‐related strain within the MDZ and serves as an example for upper crustal deformation when interpreting other hybrid tectonic systems.

Light detection and ranging (lidar) data were provided by Humboldt Redwood Company, LLC. Lidar of the study area is now available on the 3DEP LidarExplorer web interface https://apps.nationalmap.gov/lidar-explorer/ (last accessed November 2023). Post‐processed elevation data extracted from digital elevation model (DEM) by the authors and used to calculate average elevation per terrace surface along the Van Duzen River, California, are available for download at https://zenodo.org/records/10611192.

The authors acknowledge that there are no conflicts of interest recorded.

This study was conducted on land ancestral to Nongatl of the Eel River Athapaskan peoples. The authors would like to thank the Geological Society of America Graduate Student Research Grants program for partial funding of this project. The authors are grateful for the Humboldt Redwood Company’s generous access to light detection and ranging (lidar) data and property. Intellectual contributions from Melanie Michalak, Tyler Ladinsky, Jason Buck, Paul Sundberg, and Jim Falls improved this study. The authors thank Casey Loofbourrow, Michelle Robinson, J. Padgett, Heath Sawyer, Maddy Schriver, Steven Medina, Vanessa Davis and Toby Haskett for assisting with trench excavation, and Steve Tillinghast, Ian Pierce, Jesse Gremore, and Kelsey Conger for field assistance. The authors thank Christopher DuRoss for fruitful discussions on statistical treatment and Bayesian modeling. The authors would like to extend gratitude to Rich Koehler, and an anonymous reviewer for their thoughtful and thorough reviews.