Marine Paleoseismic Evidence for Seismic and Aseismic Slip Along the Hayward‐Rodgers Creek Fault System in Northern San Pablo Bay

Distinguishing between seismic and aseismic fault slip in the geologic record is difficult, yet fundamental to estimating the seismic potential of faults and the likelihood of multi‐fault ruptures. We integrated chirp sub‐bottom imaging with targeted cross‐fault coring and core analyses of sedimentary proxy data to characterize vertical deformation and slip behavior within an extensional fault bend along the Hayward‐Rodgers Creek fault system in northern San Pablo Bay. We identified and traced four key seismic horizons (R1–R4), all younger than approximately 1400 CE, that cross the fault and extend throughout the basin. A stratigraphic age model was developed using detailed down‐core radiocarbon and radioisotope dating combined with measurements of anthropogenic metal concentrations. The onset of hydraulic mining within the Sierra Nevada in 1852 CE left a clear geochemical and magnetic signature within core samples. This key time horizon was used to calculate a local reservoir correction and reduce uncertainty in radiocarbon age calibration and models. Vertical fault offset of strata younger than the most recent surface‐rupturing earthquake on the Hayward fault in 1868 CE suggest near‐surface vertical creep is occurring along the fault in northern San Pablo Bay at a rate of approximately 0.4 mm/yr. In addition, we present evidence of at least one and possibly two coseismic events associated with growth strata above horizons R1 and R2, with median event ages estimated to be 1400 CE and 1800 CE, respectively. The timing of both these events overlaps with paleoseismic events on adjacent fault sections, suggesting the possibility of multi‐fault rupture.


Tectonic Setting
The Hayward-Rodgers Creek fault system extends approximately 140 km ( Figure 1) from its intersection with the Calaveras fault in the south near the town of Alum Rock (Chaussard et al., 2015) to just north of Healdsburg (Hecker & Randolph Loar, 2018;Hecker et al., 2016). This fault system takes up substantial righ-lateral motion within the wide transform boundary between the Pacific and North American plates. The evolution of the northern Hayward and southern Rodgers Creek faults over the last 7 to 10 Ma has involved periods of extension, volcanism, and strike-slip basin development resulting in significant structural overprinting (McLaughlin et al., 2012). The most recently active strands of the Hayward-Rodgers Creek fault zone (Hecker & Randolph Loar, 2018;Watt, Ponce, et al., 2016) form a releasing bend that spans northern San Pablo Bay.
The Uniform California Earthquake Rupture Model (UCERF3, Field et al., 2013) divided the fault system into four parent sections (Figure 1) based on differences in fault zone geometry, geology, historical seismicity, and deformation style (e.g., creeping vs. locked). One of these section boundaries occurs just north of San Pablo Bay where the Hayward and Rodgers Creek faults converge (see purple bar, Figure 1). The geometry of the fault within San Pablo Bay is defined by relocated seismicity, gravity, tomographic and seismic data as vertical to steeply northeast-dipping (Langenheim et al., 2006;Lienkaemper et al., 2012;Parsons et al., 2003;Ponce, 2001;Watt, Ponce, et al., 2016). Shallow offset along the fault, as identified in chirp sub-bottom imagery in northern San Pablo Bay (Watt, Hart, & Kluesner, 2016;Watt, Ponce, et al., 2016) is predominantly down to the northeast ( Figure 2) with a majority of the off-fault deformation occurring within the hanging wall block.
Both the Rodgers Creek and Hayward faults are characterized by a combination of seismic and aseismic behavior (e.g., Chaussard et al., 2015;Lienkaemper et al., 2014). Therefore, deformation observed across the Hayward-Rodgers Creek fault in San Pablo Bay may be the result of either coseismic slip during an earthquake, postseismic slip just after an earthquake, interseismic creep, triggered slip, or a combination of the above. While the Rodgers Creek fault has not experienced a large surface-rupturing earthquake during historical times (since 1776), paleoseismic evidence suggests that the fault last ruptured during the mid-1700s (Hecker et al., 2005), with a mean recurrence interval for large earthquakes between 131 and 370 years (Budding et al., 1991;Schwartz et al., 1992). Creep has been detected along the Rodgers Creek fault north of Santa Rosa (Funning et al., 2007;Jin & Funning, 2017;Shakibay Senobari & Funning, 2019), but the fault is estimated to be ∼90% locked to the south McFarland et al., 2014). The most recent large earthquake on the Hayward fault occurred in 1868 and paleoseismic evidence suggests a mean recurrence time of 161 ± 65 years (Lienkaemper et al., 2010). The Hayward fault is one of the most well-studied creeping faults in the world, with observations from creepmeter, alignment array, repeating earthquake, GPS, and InSAR data (Bürgmann et al., 2000;Chaussard et al., 2015;Lienkaemper et al., 2014;McFarland et al., 2014;Nadeau & McEvilly, 1999;Schmidt et al., 2005;Shirzaei et al., 2013;Waldhauser & Schaff, 2021).
Since most of these land-based methods of detecting and measuring creep do not work as well in and around bodies of water, the distribution of creep along the submarine Hayward fault in San Pablo Bay is largely unknown.
There is a distinct horizontal surface creep rate gradient across the bay with an estimated 5.0 mm/yr along the northern Hayward fault just south of Point Pinole and no resolvable creep (0.3 mm/yr) along the southern Rodgers Creek fault at Wildcat Mountain . Repeating earthquakes identified along the Hayward fault beneath southern San Pablo Bay (Figure 1 and Figure S4 in Supporting Information S1) were used to suggest the fault is creeping at depth at a rate of between ∼4 mm/yr (Shakibay Senobari & Funning, 2019;Shirzaei et al., 2013) and ∼12 mm/yr (Waldhauser & Schaff, 2021). Here we present a multi-disciplinary approach to conducting on-fault marine paleoseismologyand present evidence for both vertical aseismic creep and coseismic  Hatem et al., 2021); thick purple bars, fault section boundaries (modified from Field et al., 2013); blue circles, locations of repeating earthquakes that are interpreted to indicate the presence of creep at depth (Waldhauser & Schaff, 2021); green stars, existing paleoseismic trench sites discussed in this paper. MH, Masonic Home (Lienkaemper & Borchardt, 1996); MV, Mira Vista (Lienkaemper et al., 1999); SP, San Pablo (this study); TG, Triangle G (Hecker et al., 2005); TL, Tyson's Lagoon (Lienkaemper et al., 2002). AR, Alum Rock; H, Healdsburg; L, Livermore; O, Oakland; SE, Sears Point; SF, San Francisco; SR, Santa Rosa. Creek fault (from Watt, Ponce, et al., 2016); blue line, axis of fault-bend basin in northern San Pablo Bay (from Watt, Ponce, et al., 2016); thick black lines, chirp seismic profiles highlighted in the text. (b) Chirp seismic image crossing the fault-bend basin along with vibracore locations straddling the fault. offset along the Hayward-Rodgers Creek fault within the bay. We discuss the implications of these results for both earthquake dynamics and earthquake hazards.

Hydraulic Mining
The hydraulic gold-mining process in the Sierra Nevada foothills mobilized a large volume of sediment that was transported through fluvial catchments and then deposited throughout San Pablo Bay between 1852, when hydraulic mining began, and 1884, when hydraulic mining was halted (Gilbert, 1917;Higgins et al., 2007;Jaffe et al., 2007;Krone, 1979). The gold extraction process introduced several contaminants, including elemental mercury (Hg), into the mining debris. The hydrothermal mineralization process that precipitates gold (Au) also produces magnetic iron sulfide minerals. The rapid deposition of large volumes of highly magnetic, contaminant-laden mining debris left a clear geochemical signal of the onset of the hydraulic-mining period within sediments of San Francisco and San Pablo Bays (Bouse et al., 2010;Hornberger et al., 1999). Downstream transport of mining debris from the Sierra Nevada foothills to San Pablo Bay is estimated to occur primarily during annual high flow periods associated with winter runoff (Krone, 1979), suggesting mining debris deposition began as early as 1852, and likely within one to two years thereafter. We utilized detailed down-core sampling of anthropogenic metal concentrations in cores 11 and 12B to identify stratigraphic horizons that reflect environmental changes within San Pablo Bay related to the onset of hydraulic mining activities in the Sierra Nevada in 1852 CE and subsequent industrialization around San Francisco Bay after the turn of the century (1900 CE).

Seismic Imaging
We combined ultra-high-resolution chirp sub-bottom imaging with targeted vibracoring and multi-disciplinary core analyses to characterize the timing and nature of recent deformation along the Hayward-Rodgers Creek fault within northern San Pablo Bay. Chirp subbottom imagery collected in 2014 ( Figure 2) revealed the presence of a shallow fault-bend basin in northern San Pablo Bay (see Watt, Hart, & Kluesner, 2016;Watt, Ponce, et al., 2016). Chirp data were acquired with an Edgetech 512 subbottom profiler using a 2-to 12-kHz sweep, with a 20-ms record length fired 6 times per second. Processing included subsampling at twice the sample interval, static corrections, predictive deconvolution, and trace mixing. Here we identify and trace four key reflection horizons (R1-R4) that represent acoustic impedance contrasts between stratigraphic units with different physical properties that cross the Hayward-Rodgers Creek fault within this fault-bend basin. Profiles crossing the fault zone within this basin show clear down to the east vertical separation, and the profile in Figure 2b exhibits apparent growth faulting where stratal thicknesses within the hanging wall block increase relative to the footwall. Seismic stratigraphy was mapped using Kingdom Suite and all sub-bottom depths and stratal thicknesses were estimated assuming 1,600 m/s sound velocity. Vertical resolution of the chirp imagery was approximately 5-10 cm.

Coring
A series of vibracores were collected in 2016 along and across the Hayward-Rodgers Creek fault to sample and date offset unconformities and distinct marker horizons. The combination of a dense grid of sub-bottom imagery and targeted coring allowed us to characterize both on-fault and off-fault deformation of key sub-bottom horizons across the most recently active strand of the Hayward-Rodgers Creek fault (Watt, Ponce, et al., 2016). During the 2016 coring effort, we collected 20 vibracores throughout San Pablo Bay using a Rossfelder P-5 vibracoring system ( Figure 2 Core lengths are corrected to reflect depth below the seafloor.

Table 1
Vibracore Sample Locations and Lengths that cross the fault bend basin. Cores 11 and 12B straddle the Hayward-Rodgers Creek fault, and core 15 is in the southern portion of the fault bend basin (Figure 2). Recovered core lengths ranged from 150 to 364 cm. All vibracores were cut into 150-cm sections and run through a multi-sensor core logger (MSCL) to take core photos and record gamma density, loop sensor magnetic susceptibility, P-wave velocity, and electrical resistivity. Core logger measurements were made on whole sediment cores at 1-cm intervals. Vibracores were then split, described, and sub-sampled for radiocarbon, radioisotopes, and anthropogenic metal concentrations. All core analysis data can be found at https://doi.org/10.5066/P9BDEB3K (Takesue et al., 2021). Cores were tied to seismic stratigraphy by hanging each core from the seafloor reflector and comparing core stratigraphy and gamma density variations to reflection character in the chirp imagery (see Section 6.1 for more detail).

Geochronology
Sediment age profiles reconstructed from historical bathymetry (Higgins et al., 2007;Jaffe et al., 2007) suggest the upper 2-3 meters of sediment within the north-central portion of San Pablo Bay were deposited since the mid-1850s. In order to construct a robust stratigraphic age model for the historical-and prehistoric-age strata in the cores, we integrated high-resolution chirp sub-bottom imagery with analysis of down-core trends in physical properties, detailed down-core radioisotope analysis, measured concentrations of anthropogenic metals, and radiocarbon dating. All dates represent CE, unless otherwise noted.
Chronologies for recent sediments (<100 years) were developed based on core depth profiles of Cs-137 and Pb-210 for two of the cores 11 and 12B in this study. Cesium-137 (Cs-137) and lead-210 (Pb-210) are radioisotopes commonly used to obtain chronologic information from depositional environments such as bays, lakes, and estuaries (e.g., Fuller et al., 1999). Lead-210 is a natural radioactive form of lead that is part of the uranium decay series. Unsupported (excess) Pb-210 results from decay of radon (Rn-222) in the atmosphere, and subsequent fallout results in a near continuous input. Radioactive decay of excess Pb-210 relative to radon in sediments provides a measure of sediment age as a function of depth for the last ∼110 years based on its half life of 22.3 years . Cesium-137 (half life 30.1 years) is a man-made radionuclide that provides time horizons related to the first occurrence (1953) and maximum input (1963) of atmospheric fallout from atomic weapons testing.
Concentrations of anthropogenic metals Ag, Al, As, Cd, Cu, Hg, Pb, Sb, Sn, and Zn were determined in cores 11 and 12B (Takesue et al., 2021) to identify stratigraphic horizons related to the onset of hydraulic mining activities in the Sierra Nevada in 1852 CE and subsequent industrialization around San Francisco Bay after the turn of the century (1900 CE). Subsamples for anthropogenic metals and other elements were collected approximately every 10 cm in Core 12B from 0 to 220 cm, and at higher resolution in strata bracketing the onset of hydraulic mining: every 2 cm in Core 12B and every 4 cm in Core 11. The fine fraction of sediment consisting of particles with diameter <0.063 mm was separated by dry sieving through stainless steel sieves and sent to AGAT Laboratories, a nationally accredited laboratory, to determine element contents. Sediment was decomposed with a four-acid digestion and quantified by inductively coupled plasma optical emission spectroscopy (ICP-OES) for major and some minor elements and inductively coupled plasma mass spectroscopy (ICP-MS) for minor and all trace elements. Geologic reference materials, internal standards, and consistency standards were used for calibration and quality assurance. Average differences of two sets of duplicate samples in pre-mining strata were 1% (As, Cu, Zn); 2% (Al); 3% (Pb); 6% (Cd); 15% (Sb); 17% (Sn); and 20% (Ag). Total mercury (Hg) was measured by cold vapor atomic absorption spectroscopy (CVAAS). All references herein to mercury contents refer to total mercury. The average Hg difference in two pre-mining duplicates was 0.02 mg/kg, or 23% of the natural background level. For comparison, at the onset of hydraulic mining, Hg increased by more than an order of magnitude to 0.23 mg/ kg. Percent recoveries in reference materials were within ±20% of known values for Hg and within ±15% for other elements.
Chronologies for sediments deposited prior to 1950 were derived from calibrated radiocarbon ages. Sampling for radiocarbon was conducted in multiple phases to focus subsequent sampling on intervals within each core that corresponded to diagnostic reflectors in the chirp imagery and/or distinct shifts in measured physical or geochemical properties. Bioturbation was prevalent throughout all cores; however, sampling was conducted to avoid bioturbation when possible. Radiocarbon samples from marine shell material, primarily gastropods, were sent to the National Ocean Sciences Accelerator Mass Spectrometry (NOSAMS) facility for radiocarbon measurement by accelerator mass spectrometry (AMS) 14C (Table 2). Calcareous benthic foraminifera, while present, were not of adequate abundance to be used for dating purposes. A total of 25 samples for radiocarbon dating were extracted from cores 11, 12B, and 15 in 2-cm sediment sample aliquots. Raw ages were calibrated and referenced to years before 1950 CE using CALIB 8.2 (Stuiver et al., 2021).

Core Description and Physical Properties
Sediment in vibracores 11, 12B, and 15 was predominantly mud with thin interbeds of fine sand observed at various levels within each core. Fine sand interbeds generally ranged in thickness from <1 cm up to 10 cm. Bioturbation was evident in all cores but was most prominent in the uppermost mud unit. Core depths in Table 1 were corrected to reflect depth below seafloor.
Distinct changes in core color, fine sand content, and magnetic susceptibility were apparent in all cores and help underpin the seismic-stratigraphic interpretation presented in Section 6. Note. All samples were from mollusk shell debris (100% marine carbon). Raw ages were calibrated using CALIB 8.20 (Stuiver et al., 2021) and Marine20 "global" curve (Heaton et al., 2020). Raw ages are reported in calendar years before 1950 CE. Local reservoir correction (ΔR) of 22 years was calculated in CALIB based on independent geochemical signal marking the onset of hydraulic mining in 1852 CE (see text for detailed explanation). a Depth corrected for gaps associated with areas of no recovery at the top of each core, as well as gaps related to core compaction following recovery, but prior to logging. b Raw age also include an added variance, which is a lab error calculated for 2016 based on the reproducibility of samples at the individual lab facility (0.0026 × fraction modern). by a relatively dense (1.5-1.6 g/cc), dark (5Y3/2), often finely laminated mud unit with numerous interbeds of very fine to fine sand throughout the lower two thirds of the cores. Within this lower unit as a whole, magnetic susceptibility values vary, but are generally less than 40 nanoteslas (nT). The relatively coarse-grained interbeds correspond to warmer colors on the computerized tomography (CT) images ( Figure 3). Core 15 contained a darker (5Y3/1) mud unit at the base of the core that was not visible in either core 11 or 12B. All cores transitioned from the dark (5Y3/2) unit to an upper unit of heavily bioturbated and often soupy mud (5Y4/2) with gamma densities of 1.4-1.5 g/cc and relatively high magnetic susceptibility values ranging from 40 to 50 nT ( Figure 3).

Radioisotope Profiles
Excess Pb-210 activity penetrated to depths of 42 and 19 cm in cores 11 and 12B, respectively ( Figure 4). The depth of excess Pb-210 activity indicates an approximate age of 1906 (2016-110 years). The apparent increase in excess Pb-210 between 70 and 100 cm in core 11 ( Figure 4a) is likely the result of transport of young material downcore by bioturbation, which is prevalent within the uppermost mud unit (5Y4/2, Figure 3). The first appearance of Cs-137 (in 1953) occurred at a depth of 32 cm in both cores 11 and 12B, with maxima (in 1963) occurring at depths of 18 cm in core 11 and 15 cm in core 12B (Figures 4c and 4d). Average sedimentation rates of 0.23 and 0.55 cm/yr were estimated from Pb-210 profiles in cores 11 and 12B, respectively. Maximum sedimentation rates Figure 3. Core photographs, computerized tomography (CT) scans, lithologic descriptions, magnetic susceptibility, and gamma density logs of select vibracores (see Table 1 and Takesue et al., 2021). Core section photographs for cores 11 and 12B were inadvertently cut off before at the base of each core section, leaving an apparent gap. The 1963The , 1953The , 1906, and 1852 time horizons denoted in black to the right of the lithologic description for core 12B are based on downcore measurements of Cs-137, Pb-210, and various anthropogenic metal concentrations (see Figure 5). Ages denoted in black italics represent median calibrated AMS radiocarbon ages from derived from Cs-137 profiles were 0.47 and 0.68 cm/yr for cores 11 and 12B, respectively. Notably, the counts for both Pb-210 and Cs-137 were quite low for these cores, which makes robust interpretation based on these measurements alone unwarranted.
Based on the proximity of cores 11 and 12B (within ∼300 m) and the apparent continuity of the shallowest reflectors across the fault that is evident in chirp imagery (Figure 2b), one would expect the radioisotope profiles and associated sedimentation rates to be quite similar. While the Cs-137 profiles were, in fact, comparable, the Pb-210 profiles were quite different. The very shallow penetration depth of excess Pb-210 in core 12B, which indicates an age of 1906 (2016-110 years), is suspect, in that it is equivalent to the Cs-137 maxima depth indicating an age of 1963. Possible reasons for this apparent disconnect include (a) transport of younger sediment downcore by bioturbation, (b) shortened age resolution for Pb-210 (∼60 years, e.g., Fuller et al., 1999), and (c) measurement error associated with radon. Therefore, we use radioisotope profile ages from core 11 for our stratigraphic age model (Section 5.5) and subsequent seismic stratigraphic interpretation (Section 6.1).

Anthropogenic Metals
Distinct shifts in concentrations of anthropogenic metals, particularly mercury (Hg) and antimony (Sb), in cores 11 and 12B most likely reflect environmental changes within San Pablo Bay related to hydraulic mining activities in the Sierra Nevada and subsequent industrialization around San Francisco Bay after the turn of the century ( Figure 5). Concentrations of both Hg and Sb showed marked increases over background concentrations at 106 and 131 cm in cores 11 and 12B, respectively ( Figure 5). This shift in metal concentrations is coincident with sharp increases in magnetic susceptibility (+15-30 nT) and an increase in gamma density in cores 11 and 12B ( Figure 3). Similar increases in magnetic susceptibility and gamma density were noted in core 15 at a comparable stratigraphic depth (Figure 3). We attribute simultaneous increases in anthropogenic metal concentrations, magnetic susceptibility, and gamma density to the influx of relatively coarse-grained, magnetic, contaminant-laden hydraulic mining debris into San Pablo Bay starting in 1852, or a few years thereafter.
Concentrations of Hg and Sb also showed distinct increases at 31 cm depth in core 12B ( Figure 5), with a general decrease in concentrations toward the top of the core. This trend is also mimicked by numerous other metals, including Ag, As, Cd, Cu, Pb, and Zn (see Takesue et al., 2021). The increase in combined metal concentrations starting at ∼40 cm and peaking at 31 cm in core 12B are likely indicative of industrialization after the turn of the century (Hornberger et al., 1999). More specifically, increases in gold (Ag) and lead (Pb) concentrations may be related to emissions from the nearby Selby smelter that was active between 1900 and 1970 (Ritson et al., 1999). Importantly, interpretations based on anthropogenic metal concentrations are compatible with timing derived from radioisotope profiles (horizontal dashed lines, Figure 5).

Calibrated Radiocarbon Ages
Two corrections were used to calibrate the radiocarbon measurements of marine mollusk shells to calendar ages, one for the level of atmospheric 14C (i.e., a globally averaged marine reservoir age; Heaton et al., 2020) and one for the local reservoir age. The correction for atmospheric 14C is well-constrained by tree-ring chronologies for the past 10,000 years (Stuiver & Braziunas, 1993). In contrast, the reservoir age correction has large uncertainties (hundreds of years) due primarily to the regional variability of upwelling in coastal waters and the variability of local inputs of river water containing potentially old inorganic carbon from groundwater (Ingram & Southon, 1996). For example, Ingram and Southon (1996) obtained a wide range of reservoir corrections (ΔR = −126 to 307) for three samples collected near Pt. Richmond, located approximately 13 km south of the core sites in San Pablo Bay. Due to the large uncertainties associated with reservoir corrections (ΔR) in San Pablo Bay, we calculate a new local reservoir correction using CALIB 8.2 based on an independent geochemical age constraint for the onset of hydraulic mining sedimentation (1852) in core 12B.
High-resolution (2-cm intervals) downcore geochemical sampling of Hg and Sb concentrations isolate the onset of hydraulic mining sedimentation to a depth interval between 131 and 133 cm in core 12B ( Figure 5, Table 2). This geochemical analysis provides independent evidence that the sediment at that depth interval was deposited in 1852, or within 1-2 years thereafter. A radiocarbon sample from the same core depth interval (130-132 cm of core 12B), with a raw age of 660 BP (Table 2), was iteratively calibrated using CALIB 8.2 until a ΔR (i.e., ΔR = 22 years) was identified that resulted in a median calendar age of 1852. The uncertainty associated with this ΔR value was estimated to be ±1 year based on the geochemical sampling interval (±1 cm) and the estimated sedimentation rate of this interval within core 12B (∼1 cm/yr, see Section 5.5). This reservoir correction falls within the range reported by Ingram and Southon (1996) but has a much lower uncertainty. While reservoir corrections may vary in time (downcore) and space (among cores), the timeframe represented in these cores is only ∼700 years, and the distance between cores is minimal (≤5 km). Therefore, we use this new ΔR of 22 years to calibrate all the radiocarbon ages reported in Table 2.

Age Model
Using a combination of radioisotope profiles, detailed down-core anthropogenic metal concentration measurements, and calibrated radiocarbon ages, we constructed an age model for core 12B (Figure 6) that provides a framework for assigning ages to particular reflectors in chirp seismic profiles and calculating vertical slip rates based on offsets of those reflectors across the Hayward-Rodgers Creek fault. Core 12B was used because it is in the hanging wall (downthrown side) block and represents a location most likely to contain a continuous record of sedimentation.
The resulting age versus depth profile ( Figure 6) primarily reflects interactions between sedimentation and tectonic subsidence throughout the last ∼900 years. Individual radiocarbon ages from core depths of 99, 105, 212, and 302 cm in core 12B were excluded as outliers, and all but the sample at 302 cm were associated with nearby evidence of bioturbation. Interestingly, and not surprisingly, the profile highlights an apparent increase in the linear rate of sedimentation associated with the onset of hydraulic mining in 1852. Linear sedimentation rates for the period prior to 1852 range from 0.25 to 0.31 cm/yr, depending on the inclusion or exclusion of individual suspect radiocarbon ages mentioned above. The linear sedimentation rate increases to 0.96 cm/yr after 1852. The hydraulic mining period between 1852 and 1887 has been associated with both an increase in sedimentation and widening of intertidal mudflats throughout San Pablo Bay (e.g., Higgins et al., 2007;Jaffe et al., 2007). The post-1852 sedimentation rate presented here is comparable to the 0.8 cm/yr rate calculated by Jaffe et al. (2007) based on bathymetric differencing in northern San Pablo Bay for the period between 1856 and 1887.

Integrated Seismic Stratigraphic Interpretation
Through integrated seismic stratigraphic and core analysis, we identified and traced four key reflection horizons (R1-R4), deposited in the last ∼700 years, that cross the Hayward-Rodgers Creek fault and extend throughout the fault-bend basin (Figures 7, 8, and Figure S3 in Supporting Information S1). In the following section we describe each of the four key reflection horizons, the bounding stratigraphy, evidence for growth faulting, and cross-fault correlations based on a combination of seismic stratigraphic relationships observed in chirp sub-bottom imagery, core geophysical properties, and core lithology (Table 3). Horizons R1-R4 and their respective ages and two sigma error estimates (Tables 3 and 4) are used to reconstruct the history of vertical motion across the fault in Section 6.2.
Horizon R1 (blue line, Figures 3 and 6-8) represents a distinct angular unconformity visible in sub-bottom imagery throughout the fault-bend basin in northern San Pablo Bay (Figures 7 and 8). The horizon R1 unconformity marks a sequence boundary in which the underlying reflectors are truncated by horizon R1 (Figure 8). Horizon R1 consistently separates a package of parallel to sub-parallel, mildly deformed reflectors (wavelengths of <100 m) from an overlying package of relatively undeformed parallel reflectors (Figures 7 and 8). The profile in Figure 8 also reveals the longer wavelength (km-scale) synclinal warping of strata associated with basin formation along the northeast side of the fault. There is significant growth (>0.5 m) of the overlying sediment package (unit 1a, Figure 7) across the Hayward-Rodgers Creek fault. Within cores 11 and 12B this unconformity is marked by a 5-to 10-cm thick zone of fine sand within an otherwise finely laminated mud matrix (Figure 3). This lithologic boundary corresponds to localized spikes in both gamma density and magnetic susceptibility (Figure 3). In core 15, this unconformity is represented by a very fine to fine-grained sand unit that caps a darker (5Y3/1) mud unit. In contrast to cores 11 and 12B, there is no magnetic susceptibility increase associated with this layer in core 15  Table 2) and radioisotope age horizons based on depth of excess lead-210, first appearance of cesium-137, and maximum cesium-137 depth; as well as the onset of hydraulic mining deposition in 1852 (see Figures 3-5). Colored lines represent marker horizons R1-R4 identified in chirp seismic profiles (see Figures 7 and 8). Red stars indicate interpreted earthquakes (E1? and E2) and associated two sigma modeled age range (see Figures S1 and S2 in Supporting Information S1 for Oxcal model results).
( Figure 3). In addition, there is a notable northeast-southwesttrending trough in the R1 surface that intersects the Hayward-Rodgers Creek fault at the southern end of the fault-bend basin. This trough feature is associated with abrupt truncation of strata on the northwest side of the trough and back-tilting of strata below R1 on the SE side of the trough (Figure 8). Speculation on the origins of this feature are discussed in Section 6.5. The age of the R1 unconformity is estimated in cores 11, 12B, and 15 from radiocarbon samples taken from directly above and within or below the very fine-to fine-grained sand unit. The median age of horizon R1 in cores 11, 12B, and 15 was estimated from OxCal sequence models (Lienkaemper & Ramsey, 2009;Ramsey, 2009) to be 1390, 1400, and 1 CE, respectively (Tables 3 and 4, Figure S1 in Supporting Information S1). We use the age calculated from core 12B (1400 CE) as the preferred age ( Figure S1b in Supporting Information S1), as this core location is assumed to have the most complete stratigraphic record.
Horizon R2 (green line, Figure 3) was initially identified based on core stratigraphy and associated downcore geochemical trends and was subsequently tied to chirp imagery (Figure 7, Table 3). Horizon R2 is recognized within all three cores based on the presence of a sharp lithologic contact between a lower mud unit (unit 1a, Figure 7) and an overlying very fine to fine sand unit (unit 2a, Figure 7) that corresponds with an increase in gamma density and a marked increase (∼20 nT) in magnetic susceptibility (Figure 3). This dramatic shift in lithology and magnetic susceptibility also coincides with a sharp increase in the concentration of both mercury (Hg) and antimony (Sb) (Figure 5), two metals closely linked to hydraulic mining debris (Bouse et al., 2010;Hornberger et al., 1999). Based on these combined lithologic and geochemical observations, we conclude that the base of unit 2a marks the onset of hydraulic mining debris deposition in northern San Pablo Bay in 1852 ±2 years. The reflection horizon R2 directly below unit 2a has been tied to chirp imagery (green line, Figures 7 and 8) and represents a continuous reflector that extends across the fault and throughout the fault bend basin. Stratal growth of unit 2a across the Hayward-Rodgers Creek fault suggests punctuated vertical deformation occurred across the fault prior to deposition of hydraulic mining debris in 1852 and after the deposition of the underlying mud unit 1a. We estimate the age of horizon R2 using an OxCal sequence models (Lienkaemper & Ramsey, 2009;Ramsey, 2009), based on a minimum age of 1854 directly above horizon R2 and a radiocarbon age just below the horizon, resulting in a median age of 1800 (Tables 3 and 4, Figure S2 in Supporting Information S1). The bracketing age of 1854 is used as a minimum age acknowledging that there could be a delay of up to a few years between the start of mining activities in the Sierra Nevada in 1852 and the onset of deposition within the bay.    Figure S3 in Supporting Information S1) along the western upthrown side of the fault. The poor reflection coherency above horizon R3 corresponds to the uppermost, heavily bioturbated, lower density, soupy mud unit within cores 11, 12B, and 15 ( Figure 3). The age of horizon R3 was estimated to be 1867 to 1872 based on the post-hydraulic mining sedimentation rate derived from the age model for core 12B ( Figure 6) and the depth uncertainty of the reflector (±2.7 cm).
The shallowest horizon, R4, marks the base of unit 4a (Figure 7) and denotes a seismic stratigraphic boundary above which there is no discernable discrete offset across the Hayward-Rodgers Creek fault (black line, Figures 3  and 6-8). There does appear to be a slight amount of downwarping of horizon R4 on the northeast side of the fault. This horizon is not associated with any notable changes in core lithology or measured physical properties. The age of this horizon was estimated according to the same methodology as R3 to be 1926 to 1930, suggesting vertical offset has been occurring across this portion of the fault zone up until that time and is perhaps still occurring, but offset is at a scale (<a few centimeters) undetectable with chirp imagery.

Vertical Displacement
Marker horizons R1-R3 and their respective ages and fault growth intervals throughout the seismic section shown in Figure 7 are used to reconstruct the history of vertical motion across the Hayward-Rodgers Creek fault. While marker horizons are traced on multiple seismic profiles throughout the fault bend basin, we only measure vertical throw on a single profile where we image what appears to be the most complete (thickest) sedimentary section (Figures 2b and 7). First, we measure the cumulative vertical separation (throw) of marker horizons across the fault by subtracting the upthrown block reflector depth from the downthrown block reflector depth (Table 4). The interpreted differential tectonic displacement (Table 4) reflects the net displacement for the time intervals between marker bed horizons and is calculated by taking the vertical separation calculated for each marker bed and subtracting the vertical separation measurement of the overlying marker bed or beds. Vertical slip rates are calculated using the tectonic displacement calculations and the associated preferred ages of the bounding horizons (Table 4). The uncertainty reported for the vertical slip rates in Table 4 reflects the combined uncertainty of the vertical separation measurement from chirp imagery (±2.5 cm) and the age range estimates for each marker horizon (R1-R3) described above in Section 6.1. We calculate a vertical slip rate of 1.6 ± 1.1 mm/yr for the period between 1400 and 1800 (preferred age of R1 and R2, Table 4) based on the interpreted tectonic displacement between R1 and R2 across the fault. Horizon R2, with a preferred age of 1800, has a cumulative vertical separation of 190 mm across the fault, with a ( Differential tectonic vertical displacement reflects the net displacement for the time intervals between marker bed horizons and is calculated by taking the vertical separation calculated for each marker bed and subtracting the vertical separation measurement of the overlying marker bed or beds. b Age range for horizons R1 and R2 are modeled two sigma age ranges and the age range for horizon R3 is estimated from the post-1852 sedimentation rate ( Figure 6). c Preferred age indicates the median modeled age rounded to the nearest decade for horizons R1 and R2 and the mean of the age range for horizon R3. d Inferred horizontal slip rate is calculated based on the assumption that the vertical rate represents 5%-25% of the horizontal rate for a near-vertical strike-slip fault.

Table 4 Slip Rate Estimates and Input Parameters
While the ratio of vertical dip-slip motion to horizontal strike-slip motion is unknown here, we estimate a range of horizontal slip rates based on vertical to horizontal slip ratios derived from paleoseismic observations along the southern Hayward fault and surface slip estimates derived from geodetic and geologic observations following the M w 7.1 2019 Ridgecrest earthquake in the Easter California shear zone. Measurements from the Masonic Home paleoseismic site along the southern Hayward fault (Figure 1) revealed vertical to horizontal slip ratio of 0.075 (vertical slip = 7.5% of horizontal) (MH in Figure 1; Lienkaemper & Borchardt, 1996). While this is only a single measurement along a section of the fault that exhibits reverse, rather than extensional offset, it provides a first-order constraint on the likely magnitude of coseismic offset expected for a Hayward-type event.
For comparison, systematic along-strike coseismic surface displacement measurements derived from both geodetic (Chen et al., 2020) and field-based (DuRoss et al., 2020) observations following the right-lateral strike-slip M w 7.1 2019 Ridgecrest earthquake revealed average vertical to horizontal surface slip ratios of 0.188 and 0.25, respectively. Vertical to horizontal slip ratios >25% have been observed along the extensional Cinarcik section of the North Anatolian fault (Kurt et al., 2013), however, the bend angle in this study (∼10°) is much less oblique than that along the NAF (∼23°) and would require less extension per unit length and a lower vertical to horizontal ratio. Based on these combined observations, we estimate horizontal slip rates in Table 4 based on an estimation that the measured vertical slip rate would be between 5% and 25% of the horizontal slip rate. Assuming the vertical tectonic separation in Table 4 represents the true throw, we estimate horizontal strike-slip rates of 1.6-8.2 mm/yr for the interval between horizons R3 and R4, 7.4-37.1 mm/yr for the interval between R2 and R3, and 6.5-32.6 mm/yr for the interval between R1 and R2 ( Table 4). The total vertical slip measured across the fault (820 mm) corresponds to an inferred horizontal slip rate of 5.4-27.2 mm/yr, which is compatible with the long-term geologic slip rate of 9 mm/yr for the Hayward fault (Field and 2014 Working Group on California Earthquake Probabilities, 2015).

Evidence for Creep
A number of studies have identified repeating earthquakes along the submerged trace of the Hayward fault in San Pablo Bay (see Figure S4 in Supporting Information S1; Shakibay Senobari & Funning, 2019; Shirzaei et al., 2013;Waldhauser & Schaff, 2021). These events rupture the same fault patch repeatedly and are thought to be loaded by surrounding fault creep (Nadeau & McEvilly, 1999). The presence of repeating earthquakes beneath San Pablo Bay suggests the Hayward fault is creeping at depth here. The observation of cross-fault deformation (≤60 mm) on seismic reflectors above horizon R3 that likely post-date the most recent surface-rupturing earthquake on the Hayward fault in 1868 (Figures 6 and 7) strongly suggests interseismic vertical creep is occurring at shallow depths in northern San Pablo Bay. Based on vertical displacement across horizon R3, we estimate the vertical creep rate to be 0.4 ± 0.02 mm/yr and infer a horizontal creep rate of between 1.6 and 8.2 mm/yr (Table 4). In general, estimates of vertical creep rates along strike-slip faults are limited. The limited measurements that exist for the Hayward fault suggest vertical creep rates of 0.6 ± 1 mm/yr based on paleoseismic trenching at the Masonic Home site (MH, Figure 1) along the southern Hayward fault (Lienkaemper & Borchardt, 1996) and 0.5-1.0 mm/yr based on InSAR time-series analysis along the southern Hayward fault where it intersects the Calaveras fault (Chaussard et al., 2015). These estimates of vertical creep along the Hayward fault are comparable to our estimate in San Pablo Bay and further support our interpretation that this vertical deformation is the result of aseismic creep. In addition, the horizontal creep rate inferred here is similar to the 4-12 mm/yr range estimated from repeating earthquakes along the Hayward fault beneath southern San Pablo Bay (Shakibay Senobari & Funning, 2019;Waldhauser & Schaff, 2021). Combined, these observations suggest the extent of significant (≥3 mm/yr) near surface creep on the Hayward fault likely extends >10 km north of Point Pinole.
Alternatively, the post-1868 vertical deformation could be the result of triggered slip from an earthquake on a nearby fault. Nearby moderate to large earthquakes that occurred between 1873 and 1938, the maximum and minimum ages for R3 and R4, respectively, include the M w 7.9 1906 Great San Francisco earthquake and the M w 5.8-M w 6.4 1898 Mare Island earthquake. Both static and dynamic stress changes resulting from nearby earthquakes have been considered as possible triggering mechanisms of aseismic fault slip (Bilham & Castillo, 2020;Du et al., 2003). For example, observations of a dramatic decrease in creep rate (from ∼9 mm/yr to nearly zero) along the Hayward fault following the 1989 Loma Prieta earthquake were associated with negative static stress changes induced by the earthquake (Lienkaemper et al., 1997). Calculations of Coulomb stress changes induced by the 1906 San Francisco earthquake show a marked decrease in horizontal shear stress along many Bay Area faults, including the northern Hayward fault (Harris & Simpson, 1998;Parsons, 2002), which would result in suppression rather than triggering of horizontal creep following that event. Simpson and Reasenberg (1994) calculated a positive normal stress change (unclamping) along the Hayward fault at Point Pinole (see Figure 14b in Simpson & Reasenberg, 1994) due to the 1906 earthquake that could have promoted vertical motion. However, it should be noted that the fault strike used by Simpson and Reasenberg (1994) is different than that mapped by Watt, Ponce, et al. (2016) beneath northern San Pablo Bay and therefore the resulting normal stress calculations may not be applicable to this study. While it has been postulated that the Mare Island earthquake occurred along the Hayward-Rodgers Creek fault system (Parsons et al., 2003), detailed archival accounts of the 1898 earthquake, together with a comparison of the intensity distributions for the Mare Island and nearby Mw6.0 South Napa earthquakes, argues for an earthquake source closer to Mare Island along the Franklin fault (Hough, 2014). Taking all these observations into consideration, we postulate that the vertical displacement between horizons R3 and R4 is likely the result of interseismic creep, possibly combined with triggered slip from the nearby 1898 Mare Island earthquake.

Evidence for Coseismic Offset
Vertical slip rates calculated from the vertical separation across the two oldest horizons, R1 and R2, are likely too large (>1 mm/yr) to be attributed to creep alone (Table 4) and combined with evidence for growth faulting, suggest a component of coseismic or triggered slip. Both horizons R1 and R2 are associated with overlying growth strata (Figures 7 and 9), but there are differences between the estimated amounts of vertical slip in each event and evidence for off-fault deformation and erosion associated with each horizon. While we fully acknowledge the natural (aleatoric) variability of slip behavior through time at any given fault location (e.g., Styron, 2019), we use these differences in seismic stratigraphy to try and distinguish between coseismic or triggered slip mechanisms.
The spatial extent of the R1 unconformity ( Figure 8b) and the distinct off-fault deformation visible in pre-R1 strata, in addition to the overlying growth strata, provide clear evidence that this horizon was at the seafloor during an earthquake. We argue that the formation of the R1 unconformity and the deformation of pre-R1 strata occurred as the result of rapid subsidence, strong shaking, and off-fault deformation of the shallow bay sediments during a surface-rupturing earthquake (E2, Figure 9) along this portion of the Hayward-Rodgers Creek fault. We estimate the age of this event to be between 1364 and 1462 (Table 4) and the amount of coseismic vertical offset using  Figure 7. Note that the coseismic displacements here likely include a contribution from interseismic creep. Following each growth interval, the strata on either side of the fault appear to become concordant. We assume creep is occurring at a relatively constant rate throughout the section, but at a rate that does not produce resolvable growth strata in the chirp imagery. (b) Plot of vertical displacement accumulation with time, showing the timing of interpreted earthquakes, with earthquake 1 (E1) being the youngest event. The variation in slope visible in the displacement-time curve reflects increasing throw with increasing stratigraphic depth in what we interpret to be a combination of growth intervals associated with interseismic creep punctuated by periods of coseismic displacement or triggered slip. Modified from Barnes and Pondard (2010) to represent a creeping fault.
the difference between the median modeled ages (preferred age, Table 4) of horizons R1 and R2 and assuming a background vertical creep rate of 0.4 mm/yr to be approximately 50 cm.
In the following paragraphs we discuss the possibility that E2 ruptured not only the Hayward-Rodgers Creek fault, but also a secondary northwest-southwest trending normal fault. The existence of normal faults between the Hayward and Rodgers Creek faults in northern San Pablo Bay was proposed by Parsons et al. (2003) based on tsunami modeling of the 1898 Mare Island earthquake and evidence of relative subsidence between the Hayward and Rodgers Creek faults during that timeframe from bathymetric differencing. While chirp seismic imaging did not reveal evidence of discrete dip-slip offset in the location proposed by Parsons et al. (2003), there is a buried ∼4 km-long southwest-northeast-trending trough in the R1 surface that is associated with abrupt truncation of strata on the NW side of the trough and back-tilting of strata below R1 on the SE side of the trough (Figure 8). The R1 unconformity and this trough feature are too young to have been formed by eustatic sea-level variations and there is no record of a sudden change in local water levels in the bay at this time. This trough feature could have been formed by channel erosion associated with the nearby Sonoma Creek (Figure 2), tectonic deformation (faulting), or a combination of the two. The general morphology of the trough and the overlying fill stratigraphy is consistent with erosion due to channelized flow from the creek and subsequent filling of the channel over time. However, if that was the case, the channel should be traceable in the subsurface imagery all the way to the Sonoma Creek outlet. Instead, the trough disappears >4 km from the outlet and is isolated to the fault bend basin. In contrast, the old (pre-dredge) Petaluma Creek channel system is visible in subsurface chirp imagery ( Figure  S3 in Supporting Information S1; Watt, Ponce, et al., 2016) and paleo-bathymetric surveys , connecting to its outlet into the bay.
Alternatively, the trough feature and particularly the cross-trough stratal relationships can be explained by the presence of a normal fault at depth (Figure 8), where down to the SE offset could have resulted in truncation of strata on the NW block and back-tilting of strata on the downthrown SE block. The orientation of this inferred normal fault mimics that proposed by Parsons et al. (2003). Parsons et al. (2003) used finite element modeling to demonstrate that extensional structures in northern San Pablo Bay are expected to trend southwest-northeast, perpendicular to model-predicted least-principal-stress directions. Based on these combined observations, we prefer a tectonic origin for the trough feature and postulate that E2 ruptured both the Hayward-Rodgers Creek fault and the intersecting normal fault. The formation of the trough via normal faulting likely promoted subsequent channelization of flow and channel fill deposition after faulting occurred. Unfortunately, the extensive gas within the bay sediments has not permitted direct imaging of a fault structure here or allowed for direct cross-fault correlation of strata. As a result, we are unable to directly measure dip-slip across this inferred structure. Lack of clear stratigraphic evidence for tectonic deformation in the overlying strata suggest this fault has not been active since the early 1400s, which would preclude this fault from being the source of the 1898 Mare Island earthquake and tsunami. An alternative tsunami source for the 1898 earthquake is discussed in Section 7.3.
While the evidence for R1 being an event horizon is strong, that is not the case for R2. While horizon R2 can also be mapped throughout the fault-bend basin, this horizon does not represent an angular unconformity, nor is it associated with the same widespread off-fault deformation as horizon R1. The interpretation of horizon R2 is also complicated by the fact that this horizon represents a rapid shift in sedimentation associated with the onset of hydraulic mining deposition. According to our age model (Figure 6), sedimentation rates triple above horizon R2. While the growth strata and vertical slip rate suggest horizon R2 could have formed coseismically, the observed differences between horizons R1 and R2 may reflect dramatic changes in local sedimentation and/or different fault behavior. For example, horizon R2 and the associated growth strata may have formed coseismically, but during an earthquake with less vertical slip than the R1 event. Alternatively, these growth strata formed as a result of triggered slip from rupture on a nearby fault or adjacent Hayward-Rodgers Creek fault section. Both of these scenarios account for the existence of growth strata and the lack of off-fault deformation associated with horizon R2. We estimate the age of this event (E1?, Figure 9), regardless of its slip behavior, to be between 1755 and 1845, the two sigma modeled age range for R2 (Table 4, Figure S2 in Supporting Information S1) and the amount of coseismic vertical offset using the median modeled age (1800) for horizon R2, assuming a background vertical creep rate of 0.4 mm/yr to be approximately 10 cm.
Integration of the vertical displacement estimates described above (Table 4), with identification and interpretation of growth intervals and off-fault deformation from chirp imagery (Figures 7-9) combined with existing regional geodetic and paleoseismic characterizations of slip behavior, help us to distinguish between seismic and aseismic modes of slip. In Figure 9 we present a conceptual model based on seismic stratigraphic interpretations in Figure 7 for the occurrence of coseismic displacement (CSD) and associated deposition of post-seismic growth strata (PS1-2) following two earthquakes (E1? and E2) on a creeping Hayward-Rodgers Creek fault. Following each growth interval (PS1-2), the strata on either side of the fault appear to become concordant. We assume creep is occurring at a relatively constant rate throughout the section, but at a rate that does not produce resolvable growth strata in the chirp imagery between CSD intervals. The sedimentation rates discussed in Section 5.5 (Figure 6) are consistent with the depth at which stratal concordance resumes after the post-seismic accommodation space is filled. The overall throw rate (slope) visible in the displacement-time curve (Figure 9) reflects increasing throw through time in what we interpret to be a combination of intervals associated with steady interseismic creep punctuated by periods of coseismic displacement (CSD, Figure 9) or triggered slip.

Creeping Faults and Maximum Earthquake Magnitude
Characterizing the spatial distribution of creep along faults is important for characterizing seismic hazard and understanding fault mechanics (Harris, 2017). For example, understanding how shallow fault creep influences the tectonic strain energy available to produce large earthquakes and how creep affects strong ground motions are active areas of research related to earthquake hazards. Current estimates of California's earthquake probabilities presented in the Uniform California Earthquake Rupture Forecast Version 3 (UCERF3; Field et al., 2013) consider fault creep as a factor that reduces the seismic moment available to produce large earthquakes. As a result, the maximum earthquake magnitude for creeping faults was reduced. Based on this logic, expanding the northern extent of creep along the Hayward fault, as suggested by results presented here, should result in a decrease in the expected maximum earthquake magnitude along the Hayward-Rodgers Creek fault. However, the assumption that creeping faults produce smaller earthquakes than locked faults was recently tested by Harris (2017) in a study that compared predicted and observed fault rupture areas for shallow creeping versus locked strike-slip earthquakes.
Results showed that strike-slip creeping fault earthquakes up to magnitude 6.6 produce similar fault-surface rupture areas and similar peak ground motions as locked faults. While this recent work suggests creeping faults do not necessarily produce smaller earthquakes than their locked counterparts, the authors noted that the behavior of earthquakes larger than magnitude 6.6 on creeping faults remains largely unknown due to a lack of observations (Harris, 2017).

Multi-Fault Rupture
UCERF3 (Field et al., 2013) divides the Hayward-Rodgers Creek fault system into four parent sections: the Rodgers Creek, north Hayward, south Hayward, and south Hayward extension (Figure 1). In the UCERF3 model, multi-section ruptures are permissible based on the demonstrated connectivity of the system. In evaluating the potential for cascading rupture through the extensional fault bend in northern San Pablo Bay, we must take into consideration that the path an earthquake rupture takes depends not only on fault geometry and connectivity but also on earthquake history and the slip behavior (creeping vs. locked) of the neighboring fault sections (Lienkaemper et al., 2012Lozos et al., 2015Lozos et al., , 2021Schwartz et al., 2012). Earlier work has documented the geometry and connectivity of the northern Hayward and southern Rodgers Creek faults (Lienkaemper et al., 2012;McBeck et al., 2017;Watt, Ponce, et al., 2016), and in the previous section we discussed the potential influence of creep on rupture magnitude. Here we discuss the evidence for paleoseismic events beneath northern San Pablo Bay in relation to regional paleoseismic data (green stars, Figure 1) and dynamic rupture modeling studies to provide further insight into spatial-temporal patterns of large surface-rupturing earthquakes on the partially creeping Hayward-Rodgers Creek fault zone ( Figure 10; Table 5).
Prior to this study, overlap in age estimates of paleoseismic event timing on the Rodgers Creek, northern Hayward, and southern Hayward fault sections (Hecker et al., 2005;Schwartz et al., 2014) suggested that these sections all ruptured together in a single cascading event between 1715 and 1776 or in separate, smaller events, within decades of each other. While radiocarbon uncertainties make it impossible to distinguish between a full three-section rupture and separate one-or two-section ruptures separated by decades, here we speculate on the likelihood of multi-fault rupture during the two events documented in San Pablo Bay. While the timing of E1 identified in this study does overlap with other events on adjacent fault sections, the median age estimate for this event (1800) is younger than that estimated on adjacent fault sections ( Figure 10, Table 5). In addition, the median age for E1 also conflicts with the lack of strong earthquakes along the Hayward-Rodgers Creek fault system since 1776, when a mission and presidio were built in San Francisco. Furthermore, the amount of vertical displacement associated with E1 in San Pablo Bay is quite small (∼10 cm) compared to that for E2 (∼40 cm), and there is limited to no evidence of significant off-fault deformation accompanying E1. Together, these observations provide only weak evidence that the vertical offset and growth strata associated with R2 is associated with coseismic motion during an earthquake, and, as such, argue against a large multi-fault cascading event in the mid-1700s. While there is not strong evidence here for a multi-fault Hayward-Rodgers Creek rupture during the mid-1700s, there is evidence suggesting the possibility of a cascading multi-fault rupture on this fault system between 1365 and 1461 (E2, Figure 10). The timing of the E2 documented here overlaps with the 2σ age estimates of the penultimate events along the Rodgers Creek and northern Hayward fault sections and E4 on the southern Hayward fault ( Figure 10, Table 5). In addition, the amount of coseismic vertical offset estimated for E2 in this study (0.5 m) is similar to that estimated for E4 on the southern Hayward fault (0.45 m, Lienkaemper et al., 2002). While this observation neither proves nor disproves simultaneous rupture of these fault segments, it shows that our estimates of vertical coseismic offset for E2 from seismic data are in line with those measured directly in a paleoseismic trench. Although the timing of the penultimate events on the Rodgers Creek and northern Hayward fault are poorly constrained, the amount of vertical displacement and associated off-fault deformation associated with E2 in San Pablo Bay is perhaps more in line with what would be expected for a large M7 to M7.4 multi-fault rupture. Figure 10. Plot showing the estimated timing of past surface-rupturing earthquakes from onshore (Hecker et al., 2005;Lienkaemper et al., 1999Lienkaemper et al., , 2002Schwartz et al., 1992) and offshore (this study) paleoseismic investigations along different sections of the Hayward-Rodgers Creek fault (see Table 5 for more details). See Figure 1 for locations of fault section boundaries and paleoseismic sites.  Figure 1 for locations of section boundaries).

Table 5 Summary of Dates of Hayward-Rodgers Creek Fault Post-1400 Surface Ruptures
Dynamic rupture simulations of large earthquakes on partially creeping strike-slip faults explore the influence that patches of aseismic creep have on earthquake behavior (Harris et al., 2021;Lozos et al., 2021) and provide further context for the interpretation of paleoseismic results in San Pablo Bay. Lozos et al. (2021) ran simulations using simplified partially creeping fault geometries to understand first-order behaviors, while Harris et al. (2021) incorporated more realistic Hayward-Rodgers Creek fault geometry, rock properties, and interseismic creep-rate patterns along the fault. Both studies demonstrated that the spatial distribution of aseismic creep is a critical factor controlling rupture behavior, with creep inhibiting through-going rupture in some cases but not in others. Harris et al. (2021) used estimates of the interseismic creep-rate pattern derived from geodetic and microseismicity data to inform initial stress conditions for dynamic rupture simulations (see Figure 7 in Harris et al., 2021). The southern portion of the Rodgers Creek fault from Santa Rosa to the bend in the fault at the intersection with the northern Hayward fault (Figure 1) is assumed to be predominantly locked (Jin & Funning, 2017). In the absence of creep-rate observations to the south of this bend, the creep-rate on the Hayward fault is tapered to ∼4 mm/yr near Point Pinole (Funning et al., 2007;Shakibay Senobari & Funning, 2019). In Harris et al. (2021), areas of lower initial shear stress and, thereby, lower initial stress drop were defined by fault patches with interseismic creep rates ≥3 mm/yr and separately by patches with creep rates ≥1 mm/yr. Results of the Harris et al. (2021) simulations showed that while rupture nucleating on both the Rodgers Creek and Hayward fault were capable of rupturing through San Pablo Bay in some instances, expanded areas of creep beneath San Pablo Bay and/or the inclusion of interfacial cohesion in the shallow subsurface (0-3 km depth) limited rupture propagation through the bay. For instance, simulations that nucleated on the Hayward fault using the 1 mm/yr creep rate contour to define the area of lower initial shear stress were unable to propagate at all. Simulations nucleating along the Hayward fault using the 3 mm/yr creep rate contour were only able to rupture through San Pablo Bay when cohesion was not included. Essentially, as the area of lower initial shear stress expands northward, the ability of a rupture to propagate through the connector fault decreases. Results presented here (see Section 6.2) suggest that horizontal aseismic creep rates may be as high as 2-4 mm/yr in northern San Pablo Bay. This northward expansion of significant horizontal creep may limit rupture propagation beneath San Pablo Bay, for ruptures nucleating on either the Hayward or Rodgers Creek faults. This assertion, while speculative, argues for expanded efforts to directly measure horizontal creep beneath the bay and onland to the north.

Tsunami Hazard
If the vertical deformation of the R1 horizon and associated southeast-northwest-trending trough feature described in this study was formed as the result of coseismic motion during an earthquake, combined motion on the Hayward-Rodgers Creek fault and the inferred intersecting normal fault could have displaced the overlying seawater, potentially triggering a small local tsunami. A similar scenario was proposed by Parsons et al. (2003) to explain eyewitness observations of a small tsunami and seiche following the 1898 Mare Island earthquake. Parsons et al. (2003) used hydrographic modeling to show that pure dip-slip on a normal fault (6 km length × 18 km width) with about 0.35 m average slip in northern San Pablo Bay could have produced a tsunami with maximum modeled wave height of about 0.1 m, which is quite a bit smaller than the ∼0.6-m height reported by eyewitnesses following the 1898 Mare Island earthquake (Lander et al., 1993;Toppozada et al., 1992). The explanation for the discrepancy was that perhaps a heterogeneous slip distribution, rather than the applied uniform slip distribution, might alter the maximum wave height. Interestingly, the trough feature described here is strikingly similar in geometry to the one modeled in Parsons et al. (2003). If tectonic in origin, this feature is too old though to have been the earthquake source for the 1898 Mare Island earthquake, as suggested by Parsons et al. (2003). However, dip-slip on this inferred normal fault could have produced a small tsunami like the one modeled in Parsons et al. (2003) during an earthquake around 1800, the preferred age for E2.

Conclusions
A combination of ultra-high-resolution chirp sub-bottom imaging, targeted cross-fault coring, and multi-disciplinary core analyses were used to characterize vertical deformation and slip behavior within a fault-bend basin along the Hayward-Rodgers Creek fault in northern San Pablo Bay. We identified and traced four key seismic horizons (R1-R4), all younger than approximately 1400 CE, that cross the fault and extend throughout the basin. A stratigraphic age model was developed using detailed down-core radiocarbon and radioisotope dating combined with measurements of anthropogenic metal concentrations. The onset of hydraulic mining within the Sierra Nevada in 1852 CE left a clear geochemical signature and was marked by distinct changes in physical properties and sedimentation rate within the bay sediments that was used to calculate a local reservoir correction (ΔR = 22 years) and reduce epistemic uncertainty in radiocarbon age calibration and models. The resulting seismo-stratigraphic framework was used to calculate vertical slip rates across horizons R1-R3, distinguish between recent seismic and aseismic fault behavior, and evaluate potential for multi-fault rupture. Vertical fault offset of strata younger than the most recent surface-rupturing earthquake on the Hayward fault in 1868 CE suggests near-surface creep is occurring along the fault in northern San Pablo Bay at a rate of approximately 0.4 mm/yr. In addition, we present evidence of at least one and possibly two coseismic events associated with growth strata above horizons R1 and R2, with median event ages estimated to be 1414 CE and 1800 CE, respectively. The timing of both these events overlaps with paleoseismic events on adjacent fault sections, suggesting the possibility of multi-fault rupture, although the stratigraphic evidence of large coseismic rupture is stronger for the older event.
Based on results of recent dynamic rupture modeling studies (Harris et al., 2021), the northward extension of horizontal creep (rates >1-3 mm/yr) into northern San Pablo Bay as suggested by this study may limit rupture propagation through the bay, particularly for ruptures nucleating on the Hayward fault. This on-fault marine paleoseismic investigation fills a key gap in the existing regional paleoseismic history of the Hayward-Rodgers Creek fault system (Schwartz et al., 2014), one of the most hazardous fault systems in the United States.